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EDITOR’S PREFACE 


This volume contains the papers which were presented at the Fifth Sym¬ 
posium in Applied Mathematics of the American Mathematical Society held 
at the Carnegie Institute of Technology on June 16 and 17, 1952. The sub¬ 
ject of the Symposium was Wave Motion and Vibration Theory, and the four 
sessions were devoted to Stability of Fluid Motions, Hydrodynamic Waves, 
Diffraction and Scattering Problems, and Vibration Theory. 

One of the papers appears as an abstract, because of prior arrangements for 
publication. 

All who participated in the Symposium are indebted to the McGraw-Hill 
Book Company, Inc., which, beginning with the Proceedings of the Symposium 
on Elasticity, has undertaken in these uncertain times the task of bringing the 
Proceedings of these Symposia on Applied Mathematics to the scientific public 
in book form. 

The Editor gratefully acknowledges the invaluable help afforded by the 
Committee on Arrangements, consisting of W. M. Whyburn, P. Chiarulli, P. 
Gustafson, G. H. Handelman, R. C. Meacham, L. E. Malvern, D. Moskovitz, 
and E. A. Whitman. The generous financial support provided by the Car¬ 
negie Institute of Technology greatly enhanced the success of the Symposium. 

Albert E. Heins 
Editor 
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BY 

C. C. LIN 

I. Introduction 

1. General formulation of the stability problem. The existence of two 

types of motions of a viscous fluid—laminar and turbulent—immediately 

raises the question: Which type of motion is the more likely to occur? It has 

now been generally recognized that turbulent motion is the more natural state 

of fluid motion, and laminar motion occurs only when the Reynolds number 

is so low that deviations from it are liable to be damped out. For certain 

types of flows, it has been found possible to keep the flow laminar for higher 

and higher Reynolds numbers by keeping the disturbance smaller and smaller. 

I he question then may be asked, for a given flow: Is it stable relative to 

infinitesimally small disturbances? This is the problem of hydrodynamic 
stability. 

Mathematically, the problem is as follows: Suppose the system of hydro¬ 
dynamic equations has a time-independent solution 

t 1 ) Ui{x k ) } p(x k ), T(xh) 

for the components of velocity, pressure, and temperature. Consider an 
initial-value problem with these variables slightly different from this steady- 
state solution. If the solution approaches solution (1) as time t —> oo, the 
motion is stable. Otherwise, it is unstable. It should be noted that insta¬ 
bility does not necessarily lead to turbulent motion; it could lead to another 
state of laminar motion. 

To solve the problem of hydrodynamic stability, one must therefore follow 
the solution of a system of nonlinear (quasi-linear) partial differential equa¬ 
tions, and the task is in general very difficult. In the usual approaches to the 
problem, the mathematical formulation is cast in a different way. We assume 
that, for small disturbances, the equations for the small disturbances may be 
linearized; i.e., terms quadratic or higher in the disturbances and their deriv¬ 
atives may be rejected. The resultant linear system of equations contains 
time t only through derivatives with respect to t, and hence solutions containing 
an exponential time factor e"‘ may be expected. We then have a characteristic- 
value problem, with a as the parameter. If all the characteristic values of a 
have negative real parts, the motion is stable with respect to infinitesimal 
disturbances. If some of the characteristic values a have positive real parts 
the motion is said to be unstable. The argument for this statement is briefly 

‘ST ™ te ™} WaS in f r , ted after the Presentation of this paper at the meeting. This 
added work was done with the assistance of Contract No. N5 ori-07872 of the Office of Naval 
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as follows: If the flow is given a disturbance, finite or infinitesimal, we may 
conclude that it will, in general, not die out completely. For if a disturbance 
is to die out, it must first become so small that the linearized equations are 
applicable. The small disturbance may in general be expected to contain an 
unstable mode, and it therefore cannot disappear completely. 

The validity of the process of linearization has sometimes been questioned 
in the problem of hydrodynamic stability, even though it is an often used 
process. In the literature, one finds claims of proof of complete stability in 
the nonlinear theory in cases where the linear theory definitely indicates 
instability. Careful examination of these “proofs” indicates that the claims 
are not justified. 

Even if the linearization is accepted, very complicated mathematical 
problems often have to be clarified before the characteristic-value problem can 
be solved. It is therefore no great wonder that there has been a great deal of 
disagreement among workers in this field. Some of them, unfortunately, 
result from misunderstandings. In the present article, we shall try to present 
a discussion of some of these controversial points. The two following parts 
will outline the linearized theory and the conclusions derived therefrom. This 
will be followed by a discussion of the problem of friction layers and friction 
regions. 

II. Stability of Flows Represented by Exact Solutions of 

Navier-Stokes Equations 

2. Some exact solutions and their stability characteristics. To avoid the 
complication in the various approximations involved in the theory of the 
stability of the boundary layer, we shall first consider the stability problem for 
exact solutions. There are only a limited number of such solutions in the case 
of incompressible fluid, and the following types have been repeatedly studied 
for their stability characteristics: 

(a) Poiseuille flow 

( b ) Couette flow 

(c) Plane Poiseuille flow 

(d) Plane Couette flow 

No indication of instability has ever been found for cases (a) and (d). Case (6) 
is the classical case solved by G. I. Taylor [1], who calculated its stability 
characteristics and also confirmed his theoretical results by experiments. 
Some extension of the theoretical work in this case has been done by Meksyn 
[2]. It is recognized that in this case instability does not directly lead to 
turbulent motion. 

The case of plane Poiseuille flow has been studied by many authors with 
contradictory conclusions. Heisenberg [3] was the first to reach the conclu¬ 
sion of the instability of the flow for sufficiently high Reynolds numbers, but 
he did not arrive at a critical value beyond which instability begins. The 
present writer carried out the detailed calculations [4], following the method 
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given by Heisenberg, and arrived at a Reynolds number of 5,300 based on the 
maximum velocity at the center of the channel and its half width (Fig. 1). 

Recently, Pekeris [5] again concluded that the plane Poiseuille flow is stable 
with respect to small disturbances. The present writer believes that this con¬ 
clusion arises from the fact that the formula he uses for the determination of 
characteristic values is essentially valid only in stable cases. It was decided to 
use an entirely different approach to settle this disagreement. Recently, 
L. II. Thomas [6] solved the characteristic-value problem by direct calculation 
with the help of the automatic computing machine. The results confirm those 
obtained from the Heisenberg method. 



1 

Fie;. 1. Curve of neutral stability for plane Poiseuille flow. 


3. Further discussion of the case of plane Poiseuille flow. Since the 

methods lor dealing with this classical stability problem are representative of 

the methods used in the other cases considered below, we shall present the 
theory in some detail. 

Ihe basic flow has velocity components 
( 2 ) = wi(x 2 ), u 2 = 0, a 3 = 0. 

It can be easily verified that the linearized equations for a three-dimensional 
disturbance contain variables x lf x 3 , and t only in the form of partial deriv¬ 
atives. Solutions containing the exponential factor exp i\a x Xi + a z x 3 + at] 

may therefore be expected. In fact, if the disturbance is to be finite at infinity, 
it must be spatially periodic in x x and x 3 . 

^o\\, these disturbances represent waves propagating in the direction of 
(an,0,a 3 ), which makes an angle 


t = arctan 


«3 
OC I 


(3) 
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with the Xi-axis. If we introduce a rotation of coordinates in the Xiz 3 -plane 
and take the new z*-axis in this direction of wave propagation, the com- 
ponents of velocity of the basic flow would be 

(4) ul = «!(x 2 ) cos r, u\ = 0, u* = -u^xz) sin r. 

The disturbance is now propagating in the direction of the x*-axis and inde¬ 
pendent of the variable a;*. In fact, it can easily be verified that, in this 
coordinate system, the equations governing the disturbances 
would be independent of -u* and u Thus, we have essentially a two- 
dimensional disturbance corresponding to the basic velocity components 

cos r, u* = 0, u* = 0. 

It is thus the same as the basic flow in the original problem, except that the 

velocity, and hence the Reynolds number, is reduced by the factor cos r. 

Thus,, it is only necessary to consider two-dimensional disturbances in the 

linearized theory, especially in the determination of the minimum critical 

Reynolds number. A somewhat different way of deriving this result was first 
given by Squire [7]. 

Notice that the conclusion was based on the linearized equations. It is 

well know n that three-dimensional motions are important in fully developed 

turbulence and hence may be expected to be important when the nonlinear 
terms are considered. 

The linearized differential equation for the amplitude function of two- 
dimensional disturbances is (in dimensionless form) 

(6) <£ 1V — 2 a 2 <f>" + a 4 <j) + iaR[{w — c)(0" — a 2 (p) — = 0, 

where the disturbance is given by the stream function 

t = <t>(y) exp [ia(x -cl )], 

w(y) is the basic velocity distribution in dimensionless form, 

w(2/) = 1 ~ y\ 

and R is the Reynolds number of the basic flow. The boundary planes of the 
flow are given by y = ±1, where the conditions 

(7) <K±1)=0, 0'(±1)=O 

must be satisfied. We have thus a characteristic-value problem resulting in a 
condition of the type 

(8) F(a,R,c) = 0. 

For each pair of real values a and R y there is a characteristic value c, in general 
complex. If the imaginary part c ; - of c is positive, the disturbance is unstable 
according to the linearized theory. If c, < 0, the disturbance is damped. 

If Ci = 0, the disturbance is a sustained oscillation. 
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The characteristic-value problem can be solved by numerical integration. 
However, because instability may be expected only for large values of the 
Reynolds number, the solution varies very rapidly in y and very fine steps 
must be taken. The control of errors also becomes a difficult problem. These 
difficulties have, however, been overcome by L. H. Thomas. It is estimated 
by him that the amount of work carried out would require 100 years by hand 
computation. It took 2 weeks on the high-speed IBM electronic calculator. 
As mentioned before, the results agree remarkably well with those obtained 
from the method of asymptotic solution of (6), which we shall describe pres¬ 
ently. The method of numerical calculation is at present limited to relatively 
low Reynolds numbers almost barely sufficient to produce instability for a 
limited range of values of a. Advance knowledge of the analytical results 

helps to give the correct combinations of (a,R) and to ensure that the calcu¬ 
lations would not be fruitless. 

To solve the characteristic-value problem, Heisenberg obtained four solu¬ 
tions by asymptotic methods. Two of these are given by writing 

( 9 ) 0 = <t> w (y) + • • • , 

where </> (0; is the solution of the inviscid equation 

(10) (w — — a 2 </>) — w" 0 = 0. 

Let us denote the solutions by </>i(//) and </>o(f/). The other two are given by 
the formulas of the type 

(11) 0 = ^fn(y) + —/i (y) + ■ • ^ exp [ + / \/iaR(w - c) dy], 

which may be denoted by 0 3 (j/) and My), where 0 3 shall be taken to be the 

solution which decreases exponentially as we move away from the solid 

boundary. The characteristic-value problem is also split into two parts, by 
using the symmetry property of w(ij ): ’ 

(12) (a) ^ (-1) = °’ *'(-1) = 0, 0(0) = 0, 0"(0) = 0; 

(6) 0(-1) = 0, 0'(-l) = 0, 0'(O) = 0, 0"'(0) = 0. 

Instability is found only in the latter case. In fact, if <f> is the solution of the 
inviscid equation satisfying the condition 

(!3) $'( 0 ) = o, 

then the equation defining the characteristic values is simply 
(14) fr'C" 1 ) = 03 ( I) 

*(-l) 0a( —1)‘ 

However, form (11) for 0 3 is not adequate for yielding all the results required. 
Another form of the solution, expressed in terms of Hankel functions, is 



6 


C. C. LIN 


required. This is obtained by introducing the new variable 



and expressing the solution in the form 



•A = X (0) 0?) + ex (1) (v) + • • ’ 


The procedure is suggested by considering the order of magnitude of the two 

terms ^ and iaR(w - c)*" in the neighborhood of the point y e> where w = c. 

this method actually yields four solutions at once, and the particular solution 
03 is given by 






r dr, 


f = [u>'(y c )]h. 


Heisenberg did not carry out the calculations with this new form, but he did 

suggest its use and indicated the kind of result that may be expected, including 

the correct general shape of the curve of neutral stability. The detailed 
calculations were carried out by the present author. 

4. Mathematical problems associated with the asymptotic solution of the 
rr-Sommerfeld equation. The above analysis is not free from dubious 
points, which will be discussed presently. Strictly speaking, such points can 
be settled only after a thorough investigation of the way in which the asymp¬ 
totic solutions represent approximations to the exact solutions of (6). Such 
a complete theory was developed by Wasow [8], Earlier, Tollmien [9] gave a 
satis actoiy solution in the case of real variables. Other workers were also 

able to arrive at correct answers without the benefit of the complete mathe¬ 
matical theory. 

(a) The first question that may be raised is the following: While equation (6) 
should have single-valued solutions, the asymptotic solutions (9) and (11) 
are multiple-valued, with a branch point at y = y c (which will be referred to 
as the critical point). The question arises: Which branches of these multiple¬ 
valued solutions are correct? A practical form of the answer is as follows: If 
the value of w'(y) is essentially positive at the critical point, the path for deter¬ 
mining the correct branch of these solutions should be taken below the critical 

point; if negative, above it. The rigorous form of the answer is expressed in 
terms of the real part of the integral 


Q — jj y/i(w — c ) dy. 

If one follows a path along which Rl(Q) is monotonically increasing, the 

asymptotic solutions (9) and (11) would remain valid alone; the whole path 
(see Fig. 2). 

(b) The second question is the following: Solutions of the type (9) and (11) 
are adequate for fixed values of y while aR becomes large. On the other hand, 
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the type of solution given by (16) is adequate for fixed y, as e becomes small 
{<xR becomes large). Can these two types of solutions be used at the same 
time? 

Formally, one can get some assurance that this is correct by expanding 
solutions (16) asymptotically and comparing them with solutions (9) and (11). 
But this is not strictly satisfactory from the mathematical point of view. In 
particular, the asymptotic expansion of (17) does not agree with the asymptotic 
solution (11) except in the immediate neighborhood of y c . This can actually 
be remedied [9] by redefining f in (17) by 





\/ aR(w — c) dyj 


and attaching a suitable factor to (17). Numerically, this modification turns 
out to be not significant in the present case. We may, however, note that this 
type of modification will be important for supersonic boundary layers. 


/ 

/ 

/ 



fiG. 2. Typical path T along which the asymptotic solutions remain valid. Rl(Q) = 0 
along C h Ci, C 3 ; Rl(Q ) = constant along C*, C\, C*. 


The real difficulty arises with one of the solutions of equation (10) (often 
denoted by </> 2 ). The application of the usual theory of differential equations 
to the study of the nature of the solution of (10) in the neighborhood of y = y c 
was made b} r Tollmien [10]. It is found that one of the solutions is regular, 
the other has a logarithmic singularity with derivative becoming infinite as 
V ~~* Vc- These solutions were given by Tollmien in the form 



</>, = u + ay 2 + • • • , 

w ,r 

<t>2 = 1 + by 2 + ■ • • + <t> x log y, 

c 


where a and b are constants and w' c = w'(y c ), w'J = w"(y c ). The in viscid 
solution <t > 2 obviously fails to represent a solution of (6) in the immediate 
neighborhood of y = y c , and one may expect viscous effects to be called into 
play. The point y = y c is often referred to as the “inner-friction layer.” 
This would have made no difference to the boundary-value problem, if the 
boundary points were all at finite distances (on the y-scale) from the point y c . 
Actually, however, the calculations of the neutral curve in the present case 
indicate that this distance approaches zero on the y-scale if aR —> oo. I n fact, 



8 


C. C. LIN 

JS® l0 r/' branch °!, the n u eutral curve > i1; is always at a finite distance on the 
the tvne (Q \!? ? at there ' s no guarantee that asymptotic solutions of 

asymptotic solutions (.1) „ not valid for L compute o^the twer 
branch of the neutral curve. For this reason Heisenberg suggested the use of 

L ; °u ' eV Z’ an estlmate of the error in and (which are needed in 

Li let be 0btamad . from ] the recent work of Tollmien [9] (for real y and real 
parameteis «, A, and c), and Wasow [8] (for general values of y and the param¬ 
eters). It can be shown that for „ - y c = 0(e)| that is, finite „, the per- 
centage error in <f> 2 is of the order of e and that in *' is of the order of 1/log 6 
h a small coefficient. This explains why the use of (10) gives the correct 
approximate answer. Tollmien also gave an improved solution for real c 
with percentage error of the order of e in <&. Within this error, Wasow’s 
corresponding solution for general values of c, is in agreement with Tollmien’s 

n the special case these solutions are, however, somewhat complicated to 
use for actual calculations. 

(c) There is a third problem which is not of direct interest for the calcu¬ 
lations of the neutral curve discussed above but which is of great interest 
m the general theory. What is the behavior of the solution in the limit 

“ — °? eS g ° mt0 a solutlon of the inviscid equation (10), satisfying 
inviscid boundary conditions 0(-l) = 0, *'(0) = 0? This question will be 
discussed more fully in Part IV. 


III. Stability of Flows Approximately Two-dimensional and Parallel 

5 General remarks. In the present section, we shall chiefly be concerned 
with the stability of basic flows which are parallel or nearly parallel. The 
principal types of exactly parallel flows are those mentioned in Sec. 2. How¬ 
ever, if one considers flows which are nearly parallel, there is the large class of 
flows of the boundary-layer type. Indeed, this generalization opens up a 
whole field of studies of the stability problems in a compressible fluid. Another 
striking example is the stability of the jet stream in the upper atmosphere. 

Oliver, in any such discussion it is necessary to justify the application of 

the theory of stability of parallel flows. While it has been found that this 

theory is applicable to boundary-layer flow over a flat plate and over a convex 

surface, the stability characteristics of the boundary layer over a concave 

surface and that of the pressure flow through a slightly curved channel are 

chiefly controlled by disturbances of the type found by G. I. Taylor in Couette 

flow 11]. In either case, the basic flow is, however, not much influenced by the 
curvature of the solid boundaries. 

As in the discussion in Sec. 3, the equations for small periodic three-dimen¬ 
sional disturbances can be reduced to those for two-dimensional disturbances. 
So can the boundary conditions. The transformation in the case of the 

2 y = ° 1S , far a "' a >' from point y = y c , and hence derivatives needed at that 

point are not subjected to the above complication. 
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boundary layer in a compressible fluid has been given recently by Dunn and the 
present writer [11]. However, there are two difficulties which prevent the 
neglect of the effect of three-dimensional disturbances, as is done in the case 
of plane Poiseuille flow. First, in the general case, the velocity distribution 
could change from station to station as we move downstream. Thus, although 
the transformation is valid from point to point, the history of a three-dimen¬ 
sional disturbance has no obvious correspondence with that of a two-dimen¬ 
sional disturbance. Second, in the case of the boundary layer in a compressible 
fluid, the temperature and velocity distribution in the equivalent two-dimen¬ 
sional problem does not correspond to that of an ordinary boundary layer, 
t hus, v hile the solution of the mathematical problem for a three-dimensional 
disturbance can be treated by the same method as for a two-dimensional 
one, there is no simple physical interpretation of the correspondence property. 

With the introduction of general profiles, the Rayleigh-Tollmien criterion 

for the stability of parallel flows takes on added significance. This criterion is 

as follows: If the velocity distribution does not show a point of inflection, the 

motion is stable without the influence of viscosity. On the other hand, under 

certain general restrictions, the existence of the point of inflection implies 

instability. Such a case does not exist if only the exact solutions discussed in 
Sec. 2 are considered. 

The present writer interpreted the occurrence of the point of inflection as the 

extremum of vorticity [4] and described a physical process to interpret the 

Rayleigh- lollmien criterion. Further discussions of this point will not be 
given in this paper. 

Lees and the present writer [12] generalized the condition to the case of the 

compressible fluid. In this case, the condition is that the gradient of the 

product of density times vorticity shall vanish. No simple physical interpre¬ 
tation has yet been found. 

There is one further complication in the case of nearly parallel flows The 
relative importance of the departure from parallel flow as compared with the 
viscous effect is often something to be borne in mind. For example, in the 
case of pressure flow through a channel, the calculations of Dean [13] indicate 
that a slight curvature in the channel (with a radius of curvature of 100 
channel widths, say) would cause the three-dimensional disturbance to occur 
at a lower Reynolds number than that associated with the stability of the 
straight-channel flow. In the case of flows of the boundary-layer type, the 
rate of growth of the boundary layer, and consequently the x-dependence of 
the main flow, is of the order of 1/R, where R is the Reynolds number based 
on the boundary-layer thickness, occurring in the stability equation. If such a 
dependence is neglected, as is usually the case, there may be some doubt of the 
legitimacy of including the viscous terms in equation (9). 

In the next three sections, we discuss three types of problems each with some 
distinctive characteristics. This is not intended to be a treatment of all the 
types of problems of consequence, as the previous discussions clearly indicate. 
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6. Boundary layer over a flat plate. The method of solution given originally 
by Heisenberg and described in Sec. 3 can be applied to the calculation of the 
stability characteristics of the boundary laj^er, the essential difference being a 
modification of the boundary conditions. Such calculations were performed by 
Tollmien [10] and Schlichting [14], with some improvement of the mathematical 
theory. In carrying out the calculation, Tollmien and Schlichting used solu¬ 
tions of the in viscid equation (10) in power series of y. In his problem, 
Heisenberg suggested series solutions in powers of the parameter a 2 . The 
result of Tollmien has been reproduced by the present writer by using the form 
of solution suggested by Heisenberg. The experimental work of Schubauer 
and Skramstad [15] confirms the theoretical calculations. 

Studies of the boundary layer in a compressible fluid have been carried out 
by Lees and the present writer [12]. The problems are more complicated than 
in the incompressible case. However, it is found that the boundary conditions 
on the temperature fluctuations do not enter into the characteristic-value 
problem, and the final calculations can be carried out in a manner somewhat 
similar to the incompressible case. The main difference between the two 
cases is that the disturbances can travel with a wave speed supersonic relative 
to the free stream. These disturbances are assumed to be unimportant for 
stability considerations. Subsonic disturbances, on the other hand, do not 
always exist when the free stream is at a high enough supersonic Mach number. 
In fact, Lees [16] was able to show that when heat is withdrawn from the wall 
at a sufficient rate, the boundary layer could be completely stabilized with 
respect to all subsonic disturbances. This has some experimental confirma¬ 
tion and has caused great interest and much further work. 3 

There is, however, one modification necessary to arrive at more reliable 
numerical results in all these calculations. At high free-stream speeds, the 
density distribution cannot be adequately approximated in a simple manner, 
as is done in the usual theory. This point is important especially in the case of 
high values of the wave speed c, which is the case for subsonic disturbances in a 
supersonic stream. It is then necessary to modify the solutions of the type (16) 
along the lines discussed in Sec. 4. This has been carried out by D. W. Dunn 
at the suggestion of the present writer. We record here only the solution 
which corresponds to (17) in the incompressible case. This is the one which 
will enter into the final calculations of the characteristic-value problem. The 
system of notation used is explained by Lees and Lin [12]. 

U = (ccR)-iF z (y)g 3 ^) J 

03 = (aR)-&$t(y) / </ 3 (f) df, 

J + - 

-Si = r g,(t) di, 

I 1 * 1 J + 00 

0 3 = 0 , 

3 Recent work by Dunn and the present writer indicates that there is a class of three- 
dimensional disturbances which cannot be completely stabilized by cooling. 
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where 



Fz(y) = 


$3 (y) = 

P*(y) = 

y*(y) = 
f = 



_ TTp \i(w — c) 


3m 


V\ 






w — c 

V\ 

« 2 Mif 1^3(2/), 




(y), 


fi = U 



I »)* 




w — c 
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7. Mixing of two streams and related problems. Somewhat different 

methods have to be used for the treatment of problems of “free” boundary 

layers, such as the flow in a jet, a wake, and the mixing region between two 

parallel streams. The main point is that the boundary conditions are to be 

applied at large distances from the central part of the flow field. Thus, 

asymptotic solutions of the form ( 11 ) can be used, and the boundedness of the 

solution at large distances immediately shows that they must be rejected. 

thus, viscous effects can only be found through the higher approximations of 

solution (9). This makes the nature of the stability characteristics differ 

considerably from the previous case. Studies in this direction have been 

made by Chiarulli [17], Lessen [18], Pai [19], and the present writer [20]. The 

present writer considered the mixing of two streams in the compressible case. 

Here the limitation to subsonic disturbances leads immediately to the following 

condition for stability: If the difference in the velocity of the two streams 

exceeds the sum of their velocities of sound, no subsonic two-dimensional dis¬ 
turbances can exist. 

8. Meteorological problems. One interesting application of the stability 

theory is the study of oscillations in the atmosphere. There is a strong westerly 

current, strongest at the tropopause, whose streamlines exhibit a wavy pattern. 

If one considers an idealized model of a zonal current independent of altitude, 

a two-dimensional model over a spherical earth can be studied in much the 

same way as small oscillations about an ordinary two-dimensional steady flow. 

In fact, after a Mercator’s projection, the inviscid equation for the amplitude 

of the oscillation is exactly the same as ( 10 ), except that w" is replaced by the 

gradient of vorticity, which now includes the normal component of the rotation 

ot the earth. Furthermore, since we are dealing with vorticity over a spherical 

earth, the vorticity gradient is more complicated than the second derivative 
of the velocity distribution. 

As shown by Foote and Lin [21], the stability of the flow depends on whether 

T 1 f 1 gradient van ^ s ^ es or n °L This confirms the interpretation of the Rayleigh- 

o lmien criterion given by the present writer [ 4 ] for the usual stability of two- 
dimensional parallel flows. Rossby [ 22 ] adopts the mechanism described in 

is connection for the interpretation of the movement of polar air masses. 
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Detailed discussions of the oscillations in the westerlies and their meteoro¬ 
logical implications are given by Kuo [23], 

IV. Friction Layers and Friction Regions 

9. The limiting behavior as the Reynolds number becomes infini te, We 

shall now examine in considerable detail the problem raised in Sec. 4 regarding 
the behavior of the solution of the problem of hydrodynamic stability in the 
limit of infinite Reynolds number. This question was emphasized by Heisen¬ 
berg [3], It is complicated because a formal limiting process sometimes gives 
the correct answer and sometimes does not. 

(a) Formal discussion of the limiting process. To demonstrate our point, 
we shall now carry through a formal discussion of the limiting process and indi¬ 
cate vheie such discussions lead us into difficulties. For convenience, we 
shall refer to the equations in Sec. 3, but we have in mind other cases which 
can be handled with very little modification. For example, the case of the 
boundary layer requires only the replacement of condition (13) by the condition 


*'(y) 


lim \ 
y~* « $(y) 


= —a 


If we neglect the effect of viscosity in (6), we arrive at equation (10), which 

should be solved with two boundary conditions. One of these is condition 

(13). The other is a condition at the solid wall y = - 1. Since viscosity has 
been neglected, this condition is obviously 


corresponding to u' 2 = 0. 
(14), since 


0(-l) = 0 

In a crude manner, this can be easily justified from 


>im TT7 = 0, 

aR—> * 03 V 1 ) 

whether 0 3 is given in form (11) or (17). However, if 4-' becomes infinite as 
aR —> oo, further investigation is needed. 

Supposing that the inviscid formulation can be justified (as is the case for 
the neighboihood of the first branch of the neutral curve for the boundary layer 
under an adverse pressure gradient), we have a well-defined characteristic- 
value problem, except for the fact that there is a singularity of equation (10) 
at = c (if w" ^ 0 at the same point). Thus, when viscous effects are con¬ 
sidered, they are expected to play a role at the wall and possibly at the critical 
point w = c. From the form of solutions (11) and (17), it can be surmised 
that the outer friction layer at the wall has a thickness at most of the order of 
(aR)-i, which approaches zero as the Reynolds number becomes infinite. 
Such a layer is therefore quite similar to the usual boundary layer of Prandtl. 

In the neutral case, the critical point w = c lies on the real axis, generally 
corresponding to an interior layer of the fluid. Hence, we may expect to have 
an inner friction layer in this case. In the case of amplified and damped 
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disturbances, equation (10) has no singularity on the real axis. One might 
then argue that there is no inner friction layer in these cases. In fact, equation 
(10) indicates clearly that, along the real axis, amplified and damped disturb¬ 
ances are complex conjugates. This means that to each amplified solution 
there corresponds a damped solution for the same value of a. This is difficult 
to reconcile with the first branch of the neutral curve (Fig. 3) for the boundary 
layer under an adverse pressure gradient. There are amplified solutions for 
a < a, and damped solutions lor a > <*„. By the complex-conjugate relation, 

there would be both amplified and damped disturbances for values of a larger 
and smaller than a g . How can we then ex¬ 
pect the neutral curve to separate the stable 
and unstable regions in the (a,i?)-plane? 

(b) Solution of the above dilemma. This 
difficulty disappears if we accept the follow¬ 
ing conclusions, first reached by the present 
writer in 1944. 

1. In the case ol finite amplification, the 
solution <f>(y) ol the characteristic-value 
problem tends to a solution of the inviscid 
equation, throughout the part of the real 
axis corresponding to the field of flow, as the 
Reynolds number becomes infinite. Physi¬ 
cally, this is found to mean that the effect 

of viscosity disappears completely in the interior of the fluid in the limit 

there is no inner friction layer at sufficiently large Reynolds numbers. 

2 In the case of finite damping, the solution </>(y) of the characteristic-value 

pro em does not tend to a solution of the inviscid equation along the whole 

section of the real axis of the complex ?/-plane which corresponds to the field 

ol tlow. (It does in a certain region of the complex y-plane.) The effect of 

viscosity persists at least at two layers in the interior of the fluid as the Rey- 

no ds numbers tend to infinity. The only possible exception to this statement 
is the singular case with w"(y c ) = 0. 

3. As a corollary, it is concluded that damped and amplified disturbances 

are not complex conjugates in the limiting case, as shown by an examination of 
the inviscid equation (10). 

1 hese conclusions appear rather surprising at first sight, but they are con¬ 
firmed by the work of Wasow [8]. Wasow further showed that the whole 

region between the two viscous layers in the damped case must be regarded 
as viscous. 

will ho *T CiSe ™ athematical formulation and proof of the above conclusions 

to su PP ^h“r section - we shai1 ° utiine here a simpie argument 

Wp 1 ^ inV !f Ci wu 0lUti ° n is mult 'P le - val ued (except when w"{y c ) = 0) 
We know that the correct branch is obtained by following a path Mow the 


Fig. 3. Typical form of the neutral 
curve for boundary layer under an 
adverse pressure gradient; first branch 
has an asymptotic a , 0. 
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point y c . The real axis, in the case of damped oscillations, is above y c . Thus, 
the inviscid solution cannot be a valid approximation along the whole section 
of the real axis, and viscous effect must be present. Similar arguments show 
that, in the case of amplified solutions, the inviscid solution is a valid approxi¬ 
mation along the real axis. This is essentially the original reasoning used by 
the present writer. 

In the following two sections, the limiting process aR -> oo will be discussed 
in some detail. It will be seen that there are several cases to be distinguished 
depending on the limiting behavior of the complex wave velocity c. 

10. Cases in which the critical point does not approach the end points in 
the limit. In the general study of the inviscid problem, the critical point y c > 
where w = c, does not coincide with either of the end points 2/1 and y 2 where 
boundary conditions are imposed. If y c does coincide with either 2/1 or 2 / 2 , a 
special investigation is necessary. We shall therefore first consider the asymp¬ 
totic approach of the viscous solution to those inviscid solutions for which 
y 1 “ Vc and 2/2 — 2/c remain finite while aR —> 00. 

In such cases, the characteristic-value problem can be studied by means of 
the asymptotic solutions of the types given by equations (9) and (11). It has 
been found by Wasow that the curves Rl(Q) = constant, especially Rl(Q) = 0, 
are important for the asymptotic theory. The geometry of such curves rela¬ 
tive to the real axis is shown in Fig. 4 for neutral, amplified, and damped dis¬ 
turbances (see also Fig. 2). 

In this figure, the curves C 1 , C 2 , and C 3 are the three branches of the curve 

Rl(Q) = 0 . 

The lelative position of the real axis is also* shown for all three cases. In the 

neutral case, it passes through the point y CJ where w = c. In the amplified 

case, the point y c is above the real axis by an amount I which is proportional 

to Ci for small c,-. In the damped case, the point y c is below the real axis by a 

corresponding amount. The points 2/1 and 2/2 on the real axis are the end points 

where the boundary conditions must be satisfied, with 2/1 corresponding to the 
wall. 

Clearly, other kinds of geometrical relations are possible. To fix our ideas, 
we shall limit ourselves to the kind shown in Fig. 4. The discussion of other 
possibilities will be taken up later. 

The shaded areas represent regions where the solution 0£ o) of the inviscid 
equation (10) ceases to be a good asymptotic representation of a solution of the 
complete equation (10). In fact, the solution 0 2 which has the behavior 

02 ~ <t>\ 2 0) = 1 + a(y — 2 /c ) 2 + * • • + ^~r log ( y — y c )<t> 1 

c 

in Si and S 2 diverges [8] exponentially in the sector S 3 , 


02 ~ (aR)-i exp (\/^R Q)[f 0 + 0(aR)-i] 
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where Q is defined by (18) with Rl(Q ) > 0 in $ 3 . The circular shaded region 
is associated with the fact that dfffi/dy becomes infinite as y — y 0 - 7 * 0 . 
From the solution of the form (16), it may be surmised that the radius of the 
shaded region is of the order of (aR)~t. This may be confirmed by the work 
of Wasow (cf. Sec. 4). 

These shaded areas, related to the unusual behavior of the solution <t> 2) are 
to be associated with the frictional effects in the interior of the fluid. The 



C 2 Cj 



c 2 c, 



Fig. 4. “Viscous regions" (shaded) in the 
turbances. 


C 2 Ci 



of neutral, self-excited, and damped dis- 


friction layer near the solid boundary—the outer friction layer—is given essen¬ 
tially correctly by the simple arguments in Sec. 9. 

A description of inner and outer friction layers and regions can immediately 
be obtained from Fig. 4 : 

1. The outer and inner friction layers or regions are distinctly separated 
irom each other (for otR sufficiently large). 

2. The outer layer has a thickness of the order of magnitude of (aR)~l 

approaching zero as R —> oc . y 

3. In the case of finite amplification, there is no inner viscous layer for 

sufficiently large Reynolds numbers [(a/?) - * c,]. 

4. In the case of finite damping, there is a viscous region of finite width 

however large the value of aR may become. ’ 
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, * n th . e neutral case > the thickness of viscous layer is of the order of 
(aK)-i and approaches zero as aR -> oo. 

6 . In the cases of damped and amplified disturbances with Id = _ 

and hence approaching zero with afl- =o_the disturbance has a behavior 
similar to the neutral case. 

A few comments would be proper here in view of the earlier controversy 
w ich has centered around the above conclusions. First of all, it should be 
realized that the viscous effect is to be associated with the whole third sector S, 
besides the small circular region around y c , which becomes vanishingly small as 

. ” ■ A ext ’ m the Physical interpretation, one should distinguish between 

a viscous layer whose thickness approaches zero as aR -» oc (like the boundary 

ayer of Prandtl), and a viscous region of finite width. This confusion is 
apparent in a comment by Holstein [24] on this subject. Holstein incorrectly 
stated that the present writer changed his conclusion [24] in asserting the 
existence of one viscous region. The existence of such a region does not alter 
the original conclusion of the present writer [4] that there arc two viscous layers 
in the interior of the fluid in the limit of vanishing viscosity. It corroborates 
that statement with further information about the region in between. 

It should also be noted that statement 6 does not help to answer the dilemma 

in the inyiscid case, as raised in Sec. 9. In the limit of infinite Reynolds 

number the statement concerns only the neutral case. These cases will be 

seen to have quite different characteristics. It appears that Holstein, without 

fully investigating the cases of finite damping and amplification, has used the 

results of the cases covered by statement 6 to contradict the author’s discussion 
in the cases of finite damping and amplification. 

This last point emphasizes the importance of distinguishing between different 

uniting processes Here we have the difference between finite damping and 

damping approaching zero in the manner c f = 0[(aff)H]. In the next section, 
we shall consider cases where c -* w( yi ) or w(y 2 ) as aR -> 0. 

il. Cases in which the critical point approaches one of the end points in 

e unit If the critical point y c does approach one of the end points Vl or y 2 , 
the situation is quite different. For definiteness, let us first consider the lowe^ 
and upper branches of the neutral curve in the case of the plane Poiseuille flow. 

(a) Lower branch of the neutral curve. Along this curve, the value of 


Vi = (2/1 — Vc)(aR)i 

approaches a finite limit as aR - 00 . Thus the boundary point Vl lies at a 

distance of the order of * = (aR)~i from the critical point y c , and therefore 

essentially within the circular viscous region around y c (Fig. 4). There is no 

distinction between inner and outer friction layers. There is a single viscous 

layer which extends from the solid boundary to the critical point y = y c and 
whose thickness approaches zero as aR — > oo . 

The same general conclusions hold for damped or amplified disturbances if 
the complex variable Vl approaches a finite limit as aR -» oc . It should be 
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emphasized that in the limit of infinite Reynolds number the disturbance 
approaches a steady-state deviation with c = 0. There is no distinction 
between amplified, damped, or neutral oscillations in the limit. 

However, it is possible for the characteristic value c to approach zero in a 
somewhat different manner, as happens in the next case. 

(6) Upper branch of the neutral curve. Along the upper branch, the value of 

Vi = (yi - y c )(aR)* 

becomes infinite as aR —> oo. However, y\ — y c still approaches zero. In 
fact, for plane Poiseuille motion, 

rji ~ (aR) + f* 

Vi - Vc ~ (alt.)-* 


Note that rj x is the ratio of the distance y i — y c to the thickness of the friction 
layer (radius of the circular region around the point y c ). Thus we may still 
regard the outer and inner friction layers as separated from each other. 

The above are relatively simple cases. In the case of amplified and damped 
solutions, the condition y c —> y Y implies c r —> 0 and c, —» 0. The exact manner 
in which these two limiting processes depend on the Reynolds number aR may 
still vary, and this leads to descriptions of the inner and outer friction layers 
somewhat different from the above. We shall not go into the details of such 
discussions here. 


12. Concluding remarks. The above discussions serve to illustrate the 
difficulties in the limiting process. It is necessary that these be completely 
understood; otherwise, one would wonder what happens in the limit of aR 
becoming infinite if one follows a path of constant a in Fig. 1. Obviously, the 
limiting point must correspond to a damped solution, but we know that the 
inviscid equation does not have a damped solution valid along the real axis. 
The above discussions show that the limit of a solution with finite damping 
does not satisfy the inviscid solution along the whole section of the real axis 
corresponding to the fluid field (except possibly in the exceptional case of 
w "(Vc) = 0), and consequently the effect of viscosity is not negligible in such 
limiting processes of vanishing viscosity. 

Two important reasons may be suggested for the cause of the complication 
in the discussion of the limiting process. 

1. Different limiting processes are often involved, with different results. 
This is clearly seen by contrasting the discussions of Secs. 10 and 11 and con¬ 
sidering the cases in Sec. 11. 


2. 1 he solution 0 2 has an unusual behavior, as explained in Sec. 10. With¬ 
out realizing this, one would have associated the viscous effects solely with the 
circular neighborhood of the size (aR)-* around the critical point y c . The diffi¬ 
culty in the formal limiting process, discussed in Sec. 9, would then appear. 

Finally, it may be remarked that the behavior of the solution with finite 
damping (Sec. 10) closely resembles the structure of turbulence discussed by 
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Batchelor and Townsend [25], The fluid motion exhibits a highly oscillatory 
behavior over a finite part of the field and very slow changes in another part. 

his again demonstrates an often stressed fundamental property of the motion 
of viscous fluids at high Reynolds numbers. In certain cases, the fluid behaves 
as a perfect fluid; in other cases, the effect of viscosity cannot be neglected 
even if it is very small. The increasingly refined spatial structure of the fluid 
motion is just enough to counterbalance the vanishing of viscosity and main¬ 
tain the effects of the viscous terms in the Navier-Stokes equations of motion. 

Note added in proof: 

In a recent article, Tatsumi [2G] critically discussed the method used by 

Pekeris and confirmed the work of the present writer. In two other articles 

[27], latsumi showed that the instability of the flow through a pipe could be 

caused by the instability of the laminar inlet flow prior to the formation of 
Poiseuille regime. 
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EXAMPLES OF THE INSTABILITY OF FLUID MOTION IN THE 

PRESENCE OF A MAGNETIC FIELD 

BY 


S. CHANDRASEKHAR 

1. Introduction. While the stability of two-dimensional plane flows has 
only recently been settled beyond dispute (cf. Lin [1]), there have been two 
other examples of hydrodynamic stability which though they have attracted 
much less attention have been fully understood and are probably of even 
gi eater practical significance. I he two examples are the thermal instability 
of a horizontal layer of fluid heated below and the rotational instability of 
viscous flow between rotating cylinders. In this paper we shall reconsider 
these two classical problems in the framework of hydromagnetics, i.e., when 
the fluid considered is an electrical conductor and a magnetic field is present. 
But we shall first summarize the known results in hydrodynamics. 

2. Thermal instability in hydrodynamics. The manner of the onset of 

instability in a horizontal layer of incompressible fluid has been the subject 

of investigations by Rayleigh [2], Jeffreys [3], Low [4], Pellew and Southwell 

[5], and others. The principal results established by these investigations are: 

(i) A layer of fluid heated below first becomes unstable when the Rayleigh 
number 


(i) 


It = d* 


KV 


(where g denotes the value of gravity, cl the depth of the layer considered, -/3 
is the adverse temperature gradient which is maintained, and a, «, and v 
are the coefficients of volume expansion, thermometric conductivity, and 
kinematic viscosity, respectively) exceeds a certain determinate critical value, 
(u) The motions which ensue on surpassing the critical value must have a 
cellular pattern in agreement with what had been demonstrated experimen¬ 
tally by Benard [6], The exact value of the Rayleigh number at which insta¬ 
bility sets in depends on boundary conditions such as whether a bounding 
surface is free or rigid. More particularly, the characteristic-value problem, 

the solution of which leads to the determination of the critical Rayleigh 
number, is the following: 

For a given a - (>0) find the smallest characteristic value R for which the 
equation 

(2) (£> 2 - a-yw = — a 2 RW J 'D = ~,W = W(z)] 

admits a solution which satisfies the boundary conditions 

( 3 ) W = ( D 2 - a 2 ) 2 W = 0 

19 


(for z = ±|) 
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(4) either D 2 W or DW = 0, 

depending on whether the bounding plane at z = | (similarly z = — £) is free 
or rigid. Then determine the minimum of this lowest characteristic value as 
a function of a 2 ; this minimum will specify the critical Rayleigh number R c at 
which convection through instability will first set in. 

From a direct solution of the foregoing problem it has been derived that 
(cf. Pellew and Southwell [5]) 

R c = 657.51, (both bounding surfaces free), 

(5) R c = 1,100.6, (one bounding surface free and the other rigid), 

R c = 1,707.8, (both bounding surfaces rigid). 

It may be noted here that the onset of instability by convection predicted at 
R = 1,707.8 when the layer of fluid is confined between two rigid planes has 
been confirmed experimentally by Schmidt and Milverton [7]. 

3. Rotational instability in hydrodynamics. The stability of viscous flow 
between two concentric rotating cylinders was first successfully treated both 
experimentally and theoretically by G. I. Taylor [8]. (Later investigations of 
the same problem are those of Jeffreys [3] and Meksyn [9].) In this case the 
problem is the following: The hydrodynamical equations allow the stationary 
solution 

(6) F(r) = Ar + - 

r 


for the rotational velocity at a distance r from the axis, where A and B are 
two constants related to the angular velocities fix and ft 2 with which the inner 
and the outer cylinders (of radii R\ and 7? 2 ) are rotated. Thus 



1 - mR\/R\ 
1 1 - R\/R\ 



(1 — m) 

1 1 - R 2 /Rf 


where m = fl 2 / fix. By considering a symmetrical perturbation of solution (6) 
by a disturbance characterized by a wave number X, it can be shown that the 
situation in marginal stability is governed by the equation 


(8) (DD* - X*)*» = ^ (a + v, 
where 

(9) D =4 and D* = D + 

v 7 ar r 


together with the boundary conditions 

(10) v , (DZ>* — X 2 )i/, and D*(DD* — \ 2 )v 

all vanish for r = Ri and r = /? 2 . 


INSTABILITY OF FLUID MOTION 21 

The characteristic-value problem presented by equation (8) and boundary 
conditions (10) has not been treated in its full generality. The only circum¬ 
stance under which the problem has proved amenable to treatment is when the 

difference in radii of the two cylinders is small compared with their mean i e 
when } '' 


(11) d = r 2 - Rl « i (/? 2 + Rl ) = n 0 ' 

When (11) holds, we need not distinguish between D and D* (cf. Meksyn [9]). 

And if we further restrict ourselves to the case when the two cylinders are 

rotating in the same direction, an additional simplification is possible in that 

we may replace A + (B/r») in (8) by the value O 0 of the angular velocity 

at r — R 0 . inder these circumstances the characteristic value problem 
reduces to solving 

( 12 ) (D- - a-)h = —a 2 Tv 
with the boundary conditions 


(13) 

In equation (12) 

(14) 


v = (D 2 - a 2 )v = D(D 2 - a 2 )v = 0, 


(for 2 = +£) 


T = - 


4ft„/l 


V 


2 


d 4 


is the Taylor number. It can be readily shown that this characteristic-value 

problem is identica with that appropriate for the problem of thermal instability 

in the case when both the bounding surfaces are rigid. Hence, the critical 
Taylor number for instability is [cf. equation (5)] 


(15) 


T = 1,707.8. 


A, ,vas first pointed out by Low H], this equality between the Rayleigh and 
he laylor numbers has a simple physical origin. And it may also be noted 

(,3) is ~- C* 

(16) 


du 

P ~dt P ^ U ’ S rad ) u 


mJ X H - grad p + p grad V + pvV 2 \x , 


(17) 

(18) 
and 

(19) 


curl H = 4 ttJ, curl E = - 




<9H 
dt 


div H = 0, J = <r(E +/iux H), 


dT 

dt 


+ (u • grad)? 7 = kV 2 T 


css 
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field, J the current density, p the magnetic permeability, and <r the coefficient 
of electrical conductivity. By eliminating E and J we can combine equations 
(16) to (18) to give the following pair of equations: 


(20) P 


dUi 

~dt 


and 



( 21 ) 

where 





+ pvV 2 Ui + 



dV 

dXi 




Equations (19), (20), and (21) together with div H = 0 and the equation of 
continuity are the basic equations of hydromagnetics. 

6. The problem of thermal instability in hydromagnetics. Restricting 
ourselves to the case when the impressed magnetic field ( H ) and gravity act in 
the same direction, we find that the critical Rayleigh number for the onset 
of instability depends on the strength of the magnetic field and the electrical 
conductivity, o-, through the nondimensional parameter 


(23) Q = ^ 

pv 

Indeed, we find that the characteristic-value problem which we have to solve 
for determining the critical Rayleigh number is [cf. equations (2) to (4)] 

(24) (D 2 - a 2 )[{D 2 - a 2 ) 2 - QD 2 ]W = -a 2 RW 
together with the boundary conditions 

(25) W = [(D 2 - a 2 ) 2 - QD*]W = 0 (for 2 = ±£) 

and 

(26) either D 2 W or DW = 0, 

depending on whether the bounding plane at z = ^ (similarly z = —?) is free 
or rigid. And as in the case of zero Q we must determine the minimum of the 
lowest characteristic value R of the foregoing problem as a function of a 2 . 

It has been found (Chandrasekhar [10]) that the most convenient method of 
solving the problem presented by equations (24) to (26) is by making use of the 
following variational principle: The formula 


(27) 

( +i [( DF) 2 + a*F 2 ]dz 

R - 

a 2 f** {[(Z> 2 - a*)W] 2 + Q(DTF) 2 } dz 

where 

(28) 

F = [(£> 2 - a 2 ) 2 - QD*]W, 
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provides a minimum value of It if the boundary conditions of the problem 
are satisfied and 

(29) (Z) 2 - a 2 )F = — a 2 RW 

(i.e.j if the equation governing W is satisfied). The proof of this variational 
principle will be found in Chandrasekhar [10]. But here it will suffice to note 
that it enables the solution of equations (24) to (26) by the following variational 
procedure: 

Assume for F an expression involving one or more parameters A k and which 
vanishes for z = ±£. With the chosen form of F determine W as a solution 
of the equation 

(30) [(D 2 - a 2 ) 2 - QD 2 ]W = F 

and satisfying the boundary conditions on \V at z = ±£. Then evaluate R 
according to formula (27) and minimize it with respect to the parameters A k . 
In this way we shall obtain the “best” value of R for the chosen form of F. 

The critical Rayleigh numbers listed in Table I are an extract from a more 
extensive tabulation given in [10]. 


Table I 



Rc 

Q 

Two free 

Two rigid 

One free boundary 


boundaries 

boundaries 

One rigid boundary 

0 

057.5 

1,708 

1,101 

100 

2,054 

3,757 

1,099 

500 

8,579 

10,110 

3,580 

1,000 

15,210 

17,100 

5,013 

10,000 

119,800 

124,500 

35,040 


In order to see how effective the magnetic field will be in practical cases in 
inhibiting the onset of convection, we shall consider mercury at room tempera- 
ture. We have = 1.1 X 10~ 5 , pv = 1.7 X lO' 2 , and Q = 6 X 1C 
where H is measured in gauss. For d = 1 cm. and H = 1,000 gauss, Q = 60o! 
and the critical Rayleigh numbers for the three cases considered are 

R c = 9,942, = 15.12, 

it 0 

( 31 ) Rc = 11,660, ^ = 6.80, 

it o 

Rc = 4,018, J = 3.65, 

II 0 

It would appear from these values that the predicted 
detectable in the laboratory. 


[case (i)], 
[case (ii)], 
[case (iii) ]. 
effect should be easily 
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. The problem of rotational instability in hydromagnetics. Turning next 

to the problem of the stability of viscous flow between rotating cylinders in 
J iomagnetics, we shall consider the case when the external magnetic field 

1 S "V the 2 \ dl, ' ectl0n ’ al0 »S the axis of the cylinders. In this case the 
equations of the problem (cf. Sec. 4) admit the same stationary solution (6) 

as u en no field is present. And considering symmetrical perturbations of 

this solution which are periodic in z and restricting ourselves to the case when 

<<C [C J- f^ at ion (ID] and the two cylinders are rotating in the same direc- 
ion v e find that the characteristic-value problem, the solution of which leads 

rasekha, e [i e ™ inatl ° n ^ ^ Tayl ° r nUmber for Instability, is (Chand- 


(32) 


t(^ 2 - <* 2 ) 2 + <?a 2 ]V = — a 2 T(D 2 - a 2 )^ 


together with the boundary conditions 

(33) Zty = (/)2 _ = [(£)2 _ a2)2 + = D[(Z)2 _ a2)2 + = ^ 

(for z = ±|), 

where Q has the same meaning as in equation (23). Again as in the case of 

ff i o, * we must determine the minimum of the lowest characteristic value of 
the foregoing problem as a function of a 2 . 

t has been found (Chandrasekhar [11]) that the most convenient method 
oi solving the problem presented by equations (32) and (33) is by appealing 
to the following variational principle: The formula 


(34) 


T = 


I {[(■£>* - a 2 )R] 2 + Qa 2 P 2 ) dz 


where 

(35) 


° 2 /_V ](^) 2 + a2 G 2 + a 2 Q[{DiY + aV 2 ]] dz 


P - [(D 2 - a 2 ) 2 -f- Qa 2 ]\p and G = (D 2 - a 2 )^, 


provides a minimum value of T if the boundary conditions of the problem 
are satisfied and 


(36) 


[(Z) 2 - a 2 ) 2 + Qa 2 ]P + a 2 TG = 0 


(f.e., if the equation governing \p is satisfied). The proof of this variational 

principle will be found in Chandrasekhar [11], But here it will suffice to 

note that it enables the solution of equations (32) and (33) by the following 
variational procedure: 

Assume for DP an expression involving one or more parameters Ak and which 
vanishes for z = ±£. Integrating the expression for DP, obtain P, and adjust 
the constant of integration so that P may also vanish for z = ±\. (The fact 
t at both P and DP must vanish for z = implies a certain restriction in the 
form for DP we can assume.) Next solve the equation 


(37) 


[(£>2 - a 2 ) 2 + Qa 2 ]t = p 
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for and determine the constant of integration so as to satisfy the boundary 
conditions D\p = ( D 2 — a 2 )^ = Oforz = ±^. And finally evaluate T accord¬ 
ing the formula (34), and minimize it with respect to the parameters Ak» 
In this way we shall obtain the “ 1 best” value of T for the chosen form of DP. 

The critical laylor numbers listed in Table II are an extract from a more 
extensive tabulation given in [11]. 


Table II 


Q 

T e 

Q 

T c 

0 

1,708 

1,000 

378,600 

10 

2,691 ! 

10,000 

4,459,000 

100 

17,570 j 




In order to see how effective the magnetic field will be in practical cases in 
inhibiting rotational instability, we shall again consider mercury at room 
temperature. Using the values of the physical constants appropriate for 
mercury, we find that, for d = 1 cm. and H = 1,000 gauss, Q has the value 
600; the corresponding critical Taylor number is 80 times the value for zero 
held. The predicted effect is therefore very pronounced and should be detect¬ 
able in the laboratory. Indeed, it would appear that the effect should be 
detectable even with ordinary acids and acid solutions. Thus for nitric acid 
at room temperature <r = 8 X lO” 10 , pv = 2 X 1CU 2 , and Q = 4 X 10 ~ s H 2 d 2 . 

Accordingly, under ordinary laboratory conditions the effect should become 
measurable for // lO 4 gauss. 

Finally, we may compare the inhibiting effect of a magnetic field on rota¬ 
tional and thermal instabilities. As we have already pointed out, in the 
absence of a magnetic field, the critical Rayleigh number for thermal instability 
of a layer of fluid confined between two rigid planes is the same as the critical 
Taylor number for rotational instability under the circumstances in which the 
problem has been investigated (t.e., when the cylinders are rotating in the 
same direction and d « R 0 ). When a magnetic field is present, this identity 
between the Taylor and the Rayleigh numbers no longer exists though they 
both continue to be functions of the same nondimensional parameter Q. Also 
it is seen that the inhibition of rotational instability by a magnetic field is very 
much more pronounced than the inhibition of thermal instability. 

Some related problems in hydrodynamic and hydromagnetic stability. 

1 he foregoing discussion of the classical problems of thermal instability and 
rotational instability has not included any reference to a number of related 
problems which have recently been solved. We may therefore conclude by 
briefly describing these further developments of the theory. 

(i) The problem of the thermal instability of a layer of fluid heated below 
has been extended (Chandrasekhar [12]) to include the effect of Coriolis forces 
And it has been found that when Coriolis forces are included the critical Ray- 
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leigh number for the onset of instability depends on the angular velocity of 
rotation £2 through the nondimensional parameter 


(38) 

in such a way that 



(39) R —> const. T* as T —> <x>. 

The changed dependence of 0 on y and d which this entails has bearing on a 
number of meteorological and astrophysical questions. 

(ii) The problem of the thermal instability of a layer of fluid heated below 

and subject simultaneously to Coriolis forces and an external magnetic field 

has also been investigated (Chandrasekhar [13]), and the results exhibit in a 

stiiking way the conflicting tendencies to which a conducting fluid is subject 

in its attempt to attach itself simultaneously to the vortex lines and the lines 
of magnetic force. 

(iii) the problem of the thermal instability of a fluid sphere heated within 

has been treated by Jeffreys and Bland [14] and Chandrasekhar [15], The 

analysis shov r s that the pattern of convection which is easiest to excite is a 
spherical harmonic of order 1. 

(iv) 1 he thermal instability of an incompressible sphere consisting of an 

inviscid core and a viscous mantle has also been considered (Chandrasekhar 

[16]), and it has been shown that the pattern of convection which sets in, at 

maiginal stability, in the mantle shifts to harmonics of the higher orders as the 

thickness of the mantle decreases. Thus, when the mantle extends to a depth 

of half the radius of the sphere, the harmonics of order 3 and 4 set in about 

simultaneously, while the harmonic of order 5 follows very soon afterward. 

These results have a bearing on the interpretation of the earth’s topographic 
features (cf. Vening Meinesz [17]). 

Finally it may be stated that the solution of the various problems listed 

above was possible only by suitable extensions of the variational methods 
described in Secs. 5 and 6. 
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ON FREE-SURFACE FLOWS 


P. R. GARABEDIAN 

1. The principle of minimum virtual mass. It is the purpose of this paper 

to review briefly some researches on axially symmetric cavitational flow 

which have been carried out jointly by H. Lewy, M. Schiffer, and the present 

author at the Applied Mathematics and Statistics Laboratory of Stanford 
University. 

1 he study of axially symmetric cavity flow has acquired particular interest 
with the advent of modern projectiles moving rapidly through water. The 
physical problem is characterized by motion of a projectile so rapid that the 
fluid pressure falls below the vapor pressure of water. Thus a cavity of steam 
of unknown shape forms behind the projectile. Calculation of the shape of the 
cavity is the mathematical problem to which this physical phenomenon leads. 

We study here the axially symmetric motion of an ideal incompressible fluid 

of density 1 past an axially symmetric projectile B in a cylindrical tunnel 

We assume that the motion is steady and irrotational. It can be represented 

in terms of a velocity potential <p or a stream function Let the x-axis be 

the axis of symmetry and let the xy-plane be a meridian cross section of space. 

We can represent <p and ^ as functions of x,y alone, and they satisfy the pair 
of generalized Cauchy-Riemann equations 


(i) 


Vi - 


+V 

y 


<Pt/ = — 


'tz 

y 


It is sufficient, by the above symmetry, to investigate and ^ for nonnegative 

values of y, and hence we shall consider only the half space y > 0 in the 
following. y ~ 

We require that the object B have a surface generated by rotation about the 
x-axis of a curve C which is monotonically increasing for x < 0 and monotoni- 
cally decreasing for x > 0 and is symmetric in the y-axis. We shall not dis¬ 
tinguish in our notation between a given geometrical figure in space and its 
mendian cross section in the upper half of the xy-plane. The curve C is 

assumed to intersect the x-axis at 0) and (A:„,0), and it is assumed to 

consist within the interval —k<x < k kn h >> n ^ 1 

V - h h -> n TKo , u = ’ > * > 0, Of a horizontal segment 

y - h, h > 0 I he case where C rises vertically from the x-axis is allowed 
We require that the object B lie in the tunnel 0 ^ y < h 0 , in which the flot 

1S place ’ whence it is necessary to assume that h < h 0 

e denote by L a curve, or set of curves, which bounds together with the 

egment y - h,-k g x g k, a cavity region W which lies outside B but whose 
cloM.re is contained within the tunnel 0 ^ y < h and within the strip -k < 

= ’ C den ° te by D the re 8 10n of the tunnel 0 ^ y < h which is comple- 
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mentary to B + W. It is within the region D that the above axially sym¬ 
metric flow of water takes place. B is visualized as the projectile, whereas 
W will be taken to represent a cavity of steam. 

By elimination of tp from (1), we see that \p must satisfy in the flow region 
D the second-order elliptic partial differential equation 

( 2 ) \f/ XX + yj/yy = " 

As boundary conditions on \p we require that along the tunnel wall y = ho 

h 2 

(3) * = 

while on those arcs of C and L which bound D we require 

(4) * = 0. 

On the intervals x < —k 0 , k 0 < x, of the z-axis \p is regular and vanishes. 
We suppose that \p and its derivatives are bounded at infinity. 

The virtual mass M of the axially symmetric flow \p is defined by the integral 

( 5 ) M-2,[fU 

D 

since \p — (y 2 / 2) represents the stream function of an unsteady flow in which 
B + W moves with uniform velocity against a fluid at rest at the infinite 
ends of the tunnel 0 ^ y < h 0 . We denote by 

( 6 ) V = 2tt JJ y dx dy 

b+ w 

the volume of the projectile B and cavity W. 

These preliminaries enable us to set up an extremal problem which will 
characterize the axially symmetric flow \p as a cavitational flow. We suppose 
that the projectile B is fixed, whereas the cavity W is allowed to have any 
shape lying within the rectangle —k S x ^ k ) h ^ y < h 0 . For each fixed 
fj. > 0, we seek to choose the cavity W in such a way that 

(7) M — \iV — minimum. 

The extremal problem (7) serves to determine the free-boundary curve L of W. 

Let us suppose for the moment that an analytic curve L exists bounding an 
extremal shape W for (7). We suppose that a small shift of L is made by bn 
units along the inner normal of L with respect to D. This shift generates a 
varied curve L* bounding a varied cavity W* and a varied flow region D* 
with stream function \p*. We denote by M* and V* the virtual mass of D* 
and the volume of B + W*, It is not too difficult to show [2, 10] that within 
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accuracy of the first order in bn 

(8) bV = V* — V = 2n J L y bn ds, 

(9) = M* — M = 2t J l — 1 y bn ds, 

where s is the arc length along L . But in virtue of (7) 

(10) SM — n 5V = 2tt [I (gY — 1 — n y dnds = 0 

for every choice of bn. Hence along the extremal free boundary L we must 
have 

2 

= 1 +M- 

ihe identity (11) shows that the velocity of the extremal flow is constant 
on L. On the other hand, Bernoulli’s law gives for the pressure p the relation 

(12) (^) 2 + p = const. 

Thus on L, p has the constant value 




where p 0 is the pressure at infinity. This is precisely the condition on the 

above flow necessary in order that W be a cavity of steam at constant pressure. 

Therefore our extremal flow is a cavitating flow with free boundary L and 
cavitation parameter /x. 

The earliest reference to the extremal problem (7) as a characterization of 

free-surface flow appears in the work of Riabouchinsky [10]. Existence theory 

for plane cavitational flow [2] and for axially symmetric cavitational flow [3] 

can be based on this principle. We shall proceed in the next few sections to 

indicate how a rigorous treatment of the extremal problem (7) can be carried 

through and thus how an existence theorem for axially symmetric cavitating 
flow can be obtained. 


2. Existence of a rectifiable solution. A first analysis of (7) can be made 
using the notion of symmetrization [9]. 

Symmetrization of the object B + W in the y -axis is defined to be the replace- 
ment in every meridian plane of the intersection of each horizontal line with 
B + H by a single segment of that line which is bisected by the y -axis and 
as the same length as the linear measure of the original intersection. Sym¬ 
metrization of B + W in the x-axis shall mean in this paper replacement of the 
volume of revolution B + W in space by a new volume with the property that 
its intersection with any plane perpendicular to the x-axis is a circular disk 
with area equal to the two-dimensional measure of the intersection of that 
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plane with the original volume. Symmetrization in the z-axis in space thus 
generates a certain modified symmetrization in the z-axis of the cross section 
of B + IF in the meridian xy- plane. 

Because of the symmetry and monotonicity of the boundary curve C of B, 
we see that B is unchanged by symmetrization in the z-axis or y- axis. Also, 
it is clear that the volume V of B + IF is invariant under symmetrization of 
R + B . On the other hand, it can be shown from Schwarz’s inequality or 
Minkow ski s inequality that the virtual mass M of the axially symmetric flow 
\p past B + W is either diminished or left unchanged when B + W is sym¬ 
metrized in the z-axis or in the y- axis [2, 3, 8]. 

Thus a competing shape of IF for the extremal problem (7) yields after 
symmetrization of B + IF an admissible new shape with the same value of V 
and a value of M which is either the same or smaller. Therefore it is sufficient 
to restrict the discussion to shapes such that B -f- IF is already symmetrized 
in the z-axis and in the y- axis. For a symmetrized shape of this type, the 
free curve L ascends monotonically in the second quadrant and descends 
monotonically in the first quadrant. 

1 he virtual mass M of B + IF can be estimated in terms of the virtual mass 
of a circular disk perpendicular to the z-axis with a radius equal to the maxi¬ 
mum height of L above the z-axis [3]. This estimate can be used to show that 

the quantity M — pV in (7) increases indefinitely if L rises toward the tunnel 
wall y = ho- 

We conclude that an extremal shape for IF satisfying (7) exists. For any 
symmetrized shapes II n are bounded by monotonic curves L n which satisfy 
a uniform Lipschitz condition in a coordinate system obtained by rotating 
the z-axis and y -axis through 45°. Hence, by equicontinuity, a convergent 
subsequence of shapes can be obtained from any sequence W n for which the 
associated quantities M n — pV n approach their greatest lower bound [3]. 
Ihe limiting shape IF is extremal for (7) and is bounded by a curve L which is 
monotonic in each quadrant. Furthermore, the level curves of the stream 
function \p corresponding to IF bound with the z-axis regions which are sym¬ 
metrized, and hence these level curves ascend in the second quadrant and 
descend in the first quadrant. 

3. Variational formulas. Let z 0 be any point in the interior of an arc of 
the extremal free boundary L and within the strip —k<x<k. Let p be a 
positive number so small that the circle I z — Zo | ^ 2p does not intersect B or 
the lines y = h 0 , z = —k, x = k, and let a> be a function of |z — z 0 1, with 
derivatives of all orders, which vanishes for \z — z 0 | ^ 2p and is 1 for | z — 2o| = 
p. For t in the exterior of D we set 

(13) F{z,z) = 

z — t 

For t in the interior of D , we choose rj > 0 so small that the circle \z — t\ ^ v 
is in D , and we set 
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< 14 > «M> - j4- ( . \z - <| > ,, 

< 15 > F(z,z) . J-X.1, \z — t\ < ij. 

rj~ ' 

For sufficiently small values of the complex parameter e, we consider the 
transformation 

( 16 ) z* = z + e F(z,z) 

of the independent variable z = x + ty. This transformation maps W onto 

a varied cavity W* which is admissible for the external problem (7), and it 

maps D onto a varied flow region D* with stream function We define 

further 


i**(x*,y*) = Ux,y), 

with z = x + iy and z* = x* + iy* related by (16). The functions and 

i* have the same values on the boundary of D*, and hence by Dirichlet’s 
principle 



with +** = r*(x,y), r = i*(x,y). On the other hand, by (7) 




Combining (18) and (19), we obtain 





y dx dy. 




y dx dy. 


we expand both sides of (20) in powers of e and « and use the knowledge 
that £ can be chosen arbitrarily in a small neighborhood of the origin we find 
by the usual technique of the calculus of variations 



(m + 1) 





dx dy + 4 



F thp d\p dx dy 
i dz dz y- 
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where 

l = i(±-,±\ . d \ 

dz 2 \dx dy)’ dz 2 \dx + 1 dy)‘ 

We denote by 5(0 the function which vanishes in W and is 1 in D, we denote 

by ft the intersection of D with the circle \z — z 0 1 < p, and we denote by l 

the arc of L inside the circle \z - z 0 1 < p . Letting v -» 0, we obtain from (21) 
for |< - z 0 \ < p 


(22) (/* + !) 



y dz 

z — t 



+ 4 ^ 


dx dy 


r’(z - 0 


W ««) + Q(t) = 0 , 


where Q(t ) is an analytic function in the circle |< - z 0 \ < p . Successive appli¬ 
cations of Green’s theorem to the integral over ft in (22) show in turn that 



(m + 1) 



+ 8(0 



where q(t) is continuous for \z — /| < p . 

In order to find the boundary values of \p t on l from (23), 
the function 


47r 


{t 


* 


i) 3 

+ q(l) = 0 , 


we introduce on l 



defined by integration along l. The derivatives 



exist almost everywhere on the monotonic arc l. We let z l = xi + iy x be a 
point of l where (24) holds, and we set t x = z x + f e ie , < 2 = z x - fr ie for small 
? > 0, with 6 = tt/ 4 if x x > 0 and 6 = 37r/4 if x x < 0. Thus h lies in ft, and 
/ 2 lies in IT. From (23) and an integration by parts we obtain 


(25) lim 


167T 


K 


l\—*z\ (1 +n)(t 1 - ti) 

y dz 


= lim 
^->o 



i z - 1 1 


= 2iriyii\ + lim 2 



V dz 
i z - t 2 

w(z) - 


w(z i) 


Z — Zl 


V\Z\ 


(ti - t 2 ){z - Zl ) 2 dz 

(z - 1 1 ) 2 (z - t 2 y 


One can check that the last limit appearing in (25) is zero by introducing a 
new coordinate system a, r with origin at Z\ and with the r-axis inclined at an 
angle 0 with the x-axis. Indeed, because of the monotonicity of l, the last 
integral (25) can be estimated easily [3] in terms of integrals with respect to 
a over the intervals (- oo,- € J), (-e*,-*), (-e,e), (e,c*), and ( € *,oe). 
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(26) 


lim 


y 


= —1(1 -h 


as 2 approaches L along a line inclined at 45° with the horizontal. The same 
calculation shows, with slight modifications, that i J/ z is bounded in the neighbor¬ 
hood of every point of L, except only for the intersection of L with the ?/-axis, 
near which L is not monotonic. However, the Phragm<§n-Lindelof theorem 
can be used to show the boundedness of even at that intersection [3]. The 
boundedness of the gradient of ^ together with the limit relation (26) almost 
everywhere along L constitutes a suitable generalized boundary condition 

from which we can study the free curve L without any further reference to 
the extremal problem ( 7 ). 

4. The functional equation. We turn next to a proof that the free curve L 

is analytic. Our procedure [3] is the outgrowth of earlier work on free-surface 
flows in a gravity field [ 6 ]. 

The elliptic partial differential equation ( 2 ) has a fundamental solution S 
of the form 



S(z,z;f,.0 = - 4 ( 2 , 2 ;f,f) log (z - f)(z - f) + R( Zj 2 ;f,f) 


with a logarithmic pole at 2 — f, where A and B are regular analytic functions 
of 2,2 and f,f and 


( 28 ) A ( 2 , 252 , 2 ) = 

1 he function A( 2 , 2 ;f,f) is merely the Riemann function of the partial differ- 
ential equation of hyperbolic type 


<29) *■' + 2 <r=T> *■ - 2<rb-, *■ - »■ 

which is equivalent to (2) in the real domain. We have for A the explicit 
representation [ 1 ] 


(30) 


^(M;f,f) = — z) * 


2 (f ~ f) 


F 


1 K (2-f)(2- 


2’ 2’ 1 ’ (2 - f)(z 



in terms of the hypergeometric series F. 

If T denotes a closed curve in D, we find by Green’s theorem 



2 8(fWf ,f) 

r - ? 



f) 


S(z,z ;f,f) 


d[f/ ds 
dn\ y’ 


where 5(f) now denotes the function which is 1 inside T and 0 outside T. 
By performing a suitable limit process on T, we can extend (31) to the case 
where one arc of r coincides with an arc l of L. On l we can use the formula 
( 2 b) to determine d^/dn, by the Lebesgue convergence theorem. 
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Let us write f - £ + in, T* = £ - A 7 , and let us continue the identity (31) 
analytically into the domain of two complex variables £ and i). Such an 
analytic continuation shows that 


ds 
V 

for independent values of f and f *, provided f and f* lie either both inside or 

both outside T. This identity even holds in the limiting case where [* lies 

on l and f lies either in r or outside r. We shall deal only with this latter 
situation of f*. 

The determination of the branch of log (z - f)(z - f*) which must be used 

in (32) changes by 2jrt on the arc of l between f * and f when f crosses l, but it 

does not change on the remainder of r. Thus by calculating the jump of 

each side of (32) as f crosses l and by taking into account the conditions (4) 

and (26) satisfied by \p on /, we find that for f on l, yp has continuous boundary 
values given by 




2*(r,n 

t - r* 


(1 + m)* [ A(z,z£,£*) ds. 


Also, we can differentiate (32) with respect to f and again calculate the jump 
across l Using the monotonicity of l, we can thus show that almost every¬ 
where along L the boundary condition 


(34) = (1 + + (1 + M )* J*±A(z,zt,D ds 

is fulfilled. The formulas (33), (34) merely represent the solution i/' of the 

Cauchy problem for (29) in terms of the Riemann function A, with initial 
data given along L. 

Motivated by (34), we set up for any fixed f*, with f* on L, the system of 
two integral equations 


(35) 

(36) 



_ 2_. * vKf,f*) _ ^ 

(1 + M)^(f,ff(f);f,f*) <9f f - f* J r Atf,gQ;)t,n ) 9 

</(f) = r + /_ r /(«)*&, 


for the determination, in Z) and on L near f*, of the two analytic functions 

/(f) and <?(f). ihe method of successive approximations applies to the system 

(35), (36) and can be used to prove [3, 7] the existence of a unique solution 

fyQ near f * which is analytic and is independent of the path from f * to f. All 

that is required for this existence proof is the known boundedness of the 
derivative 

We notice by (34) that along L, 


/(f) = f, 9(f) = f 
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is a solution of (35) and (36). As an integral, g(f) must be continuous in 
D + L near £*, and by the uniqueness of the solution of (35) and (36) we con¬ 
clude that the analytic function g(£) in D has on the free boundary L the values 

( 3? ) ff(f) = f- 


Let t - /(f), f = f(/) be a conformal mapping of D onto the upper half of 
the /-plane, and let 

m = m + g(m), 

*(0 = m - g(Ut)). 

Then on the interval of the real axis corresponding to L 



$(/) = real, \k(/) = imaginary, 


by (37), and hence Schwarz’s reflection principle allows us to continue $ and 
* across the real axis into the lower half plane. In particular, 



m - 


is analytic for real values of / corresponding to L, and hence the free boundary 
L is an analytic arc. 

It is now easy [3] to conclude that i is analytic on L and that the constant- 
pressure condition (11) is satisfied identically there. 

5. Remarks on uniqueness and construction. Let Zh and D 2 be two differ¬ 
ent cavitational flow regions obtained by solving the extremal problem (7) with 
two values M i and of the cavitation parameter. Suppose that we can mag¬ 
nify Dj by a factor a ^ 1 until its lower boundary curve lies just inside D 2 . 
Ihen we can show by the generalized maximum principle for (2) that 




< fa(x,y) 


in the intersection of D 2 and the magnified region D h since by our lemma on 
symmetrization the level curve 



lies above the horizontal line y = h 0 . Hence at a point of tangency of the 
fiee boundaries of Do and the magnified image of D\ 



Thus we obtain for ^ and /i 2 the useful inequality [5] 


(43) 


Ml < M2. 


From (43) we can conclude that the extremal cavitating flow obtained from 
(7) for a given value of a is unique [2, 3, 4, 5], since otherwise we would 
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arrive at the absurd inequality M It is also possible to deduce from (43) 

that the cavity T7 depends continuously on and increases steadily as /x 
inci eases [2, 3]. In particular, for one special value of the cavitation parameter 
the fiee boundary L is an analytic arc joining the ends of the segment y = h, 

k ^ x ^ k. d his statement indicates that we can prescribe the point of 
separation of the free boundary by suitable choice of /x. 


A case of special interest for the existence theorem on cavity flow completed 
above is that in which the curve C consists outside the strip —k < x < k of 
two vertical segments x = -ft and x = k, 0 ^ y g h. This is the case of a 
nose and tail consisting of circular disks. We take the point of separation of 
L fiom C to be k + ih and we set f * = ft + ih in (35) and (36). We check 
easily that for f on the segment x = ft, 0 ^ y ^ h, the expression 



9 


g *(f,n 

^ f - r 


d 



A(z,zXJ*) 


d\p ds 
dn y 


is pure imaginary. Furthermore, for pure imaginary values of g — f*, 

95) 4(r,?;r,n = real, 

( 46 ) (z,g£,S*) = imaginary 

on that segment. Hence there is a solution f,g of (35) and (36) along the 
segment x = ft, 0 ^ y g h, with 


9 — t * = imaginary, / = imaginary 

there. By the uniqueness of the solution of (35), (36) we deduce that the 
original analytic function g({) satisfying (37) has on that segment a pure 
imaginary variation. 

This information together with (38) would enable us to map the 4>-plane 
upon the ^-plane and thus determine L explicitly if we could make similar 
deductions along the x-axis. As it is, we see that on C and L the variation of 
is imaginary, while on C the variation of 4> is imaginary and on L the varia¬ 
tion of $ is real. Thus by studying the conformal mapping of the <£-plane 
onto the Sk-plane we can obtain an asymptotic development for Sk as a func¬ 
tion of 4> near the point of separation f *. This development can be used to 
show that the local behavior of the free boundary at the point of separation 
from a vertical fixed boundary is the same for axially symmetric flow as for 
plane flow, since the functions g , <£, and \k, defined by (37), have the same 
behavior in the plane case as that described above. 

Finally, the equations (33), (35), (36), and (37) allow us to construct 
explicitly an example of axially symmetric cavitational flow with an arbitrary 
given analytic curve L as free boundary. Indeed, given an analytic curve L 
we can write down explicitly an expression for the analytic function w = g(£) 
with the property (37) on that curve. For points f and f * on L we set 




(f-f*)(l + m)> 
2 


A(z y g(z)X,t*)g'(z)i dz , 

r 


(48) 
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a step which is motivated by (33), (35), and (36). We note that g({*) = f*, 
and we let f = G(w) be the inverse of w = g(£). From (48) we have 

(49) lKf,w) = —- ^ f A(z,g(z)X,w)g'(z)l dz, 

Z J G(w) 

replacing f * by w. Formula (49) is an analytic identity for f and w on L, and 
hence by analytic continuation it should hold in some regions of the z-plane 
and the tc-plane. In particular, we find [3] that 

(50) ^(f,f) = —— ^4 -- — ■ ^ [ _ Mz,g(z)£,t)g'(z)* dz 

Z J GU) 

is a real solution of (2) which represents the stream function of an actual 
axially symmetric cavitational flow with the given analytic curve L as its free 
boundary. 
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REVIEW OF SIGNIFICANT OBSERVATIONS 
ON THE MACH REFLECTION OF SHOCK WAVES 


BY 


WALKER BLEAKNEY 

The problem of the interaction of one shock front with another has received 
considerable attention in recent years and in the restricted form of the oblique 
reflection of a plane wave from a rigid wall has been the object of both theo¬ 
retical and experimental attack. In spite of the careful investigations which 
have thus far been made many features of the phenomena associated with the 
problem are not understood. The situation has been summarized in several 
recent papers [3, 4, 16]. It is the purpose of this communication to discuss the 
salient observational features of Mach reflection which must be taken into 
account in any complete theory of the phenomenon. In writing this short 
account some known points are reviewed and some new ones added, and it 
will be assumed that the reader is familiar with the main ideas and notation 
given in recent discussions [3, 16], 

Mathematical difficulties have excluded, as yet, a solution of the hydro- 
dynamical problem of a plane shock of finite strength in an inviscid but com¬ 
pressible fluid striking symmetrically a wedge of finite angle. In the absence 
of a complete solution a local three-shock theory was suggested by von Neu¬ 
mann [2] to apply in the neighborhood of the intersection of the incident, 
reflected, and Mach shock fronts. Essentially the same ideas give a solution 
for regular reflection (two-shock theory) which corresponds with the experi¬ 
mental observations. An unexpected difficulty appeared, however, when 
Smith’s results [6] indicated that the three-shock configurations were found in 
angular domains where no solutions of the three-shock equations exist. 

Stimulated by this discrepancy, a number of asymptotic solutions of the 
hydrodynamical problem valid under limiting conditions have been found. 
The work of Bargmann [9], Lighthill [19], and Ting and Ludloff [20] covers the 
case of all shock strengths at nearly glancing incidence and Keller and Blank 
[24] the case of vanishing shock strength at all angles of incidence, while 
Lighthill [26] has also treated the case of nearly head-on reflection. In so far 
as these solutions can be compared with experiment no disagreement has been 
found. Where Mach reflection appears in these asymptotic solutions, the 
reflected wave is vanishingly weak and there is no contradiction with the local 
three-shock theory. It is also worth noting that the assumption of a perfect 
fluid in the asymptotic solutions seems to be quite justified. 

The difficulty in explaining Mach reflection at finite angles has led to many 
speculations regarding the validity of the experimental measurements and 
the significance of certain observed features. So far no departure from the 
pseudo-stationary state has been detected (i.e., the whole phenomenon may 
be described in terms of Xl /t and x 2 /t, for example, where x is a coordinate and 
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good precision, and these alone suffice to show that the local three-shock theory 
does not fit. 

Most of the experimental observations have been made in air, which is 
assumed to obey the perfect-gas law. An attempt was made to see whether 
the nature of the gas had any marked influence on Mach reflection by repeating 
the experiment in argon and carbon dioxide. The first result had all the 
appearances observed in air, while the second showed large relaxation effects 
(■^S* f)i but in neither case could the results be fitted to the simple theory. 
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The treatment of the problem neglects the viscosity and heat conductivity of 
the gas, and these points obviously require examination. Unfortunately it is 
difficult expeiimentally to obtain gases with sufficient variation in these respects 
to throw any light on the question. The strongest arguments in justification 
of the assumptions of a perfect fluid are first the pseudo-stationary property 
and second the success of the two-shock theory of regular reflection without the 
necessity of introducing transport phenomena. 

Interferometer measurements of the density [16] have failed to reveal any 
angular variations such as Prandtl-Meyer expansions which might alter the 
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No clue to the puzzle in the neighborhood of the shock intersection having 
been found, it is appropriate to look at other regions of the disturbance. A 
question that arises concerns the nature of the corner where the reflection first 
begins. In some of the early work this corner was represented by the junction 
of an inclined plate with the wall of the shock tube, and investigation showed 
that there was an interaction here with the boundary layer at the later times, 
lo eliminate this disturbance as far as possible, a symmetrical wedge was used 

in all the later work, as shown in Fig. 3, but no detectable change in the reflec¬ 
tion process itself was observed. 



f 10 ' Secondary .shock in Mach reflection of very strong waves. The slip stream is 

turbulent near the wall. Note the second slip stream. Density contours are numbered. 

The densities m front of and behind the incident shock are denoted by Pl and P „ respectively. 
Angle of incidence a = 42°. J 


I he boundary conditions on the slip stream as it meets the wall are such 
that either it is tangent to the wall or its intensity vanishes there. Some 
adjustment of the stream was therefore expected, but the manner in which this 
takes place was not foreseen. As shown in the interferogram in Fig 4 the 
end of the slip stream near the wall curls up on itself in a peculiar manner, the 
effect showing up best for strong shocks of pressure ratios in the neighborhood 
oi 10. In this range the reflected and Mach shocks are quite straight near the 
triple point, and the three-shock theory fits pretty well. As the experiments 

, i. . j a point of inflec¬ 

tion, oi kink, begins to appear in the reflected wave [23] and a compression 

forms along a front extending roughly from the kink to the point at which the 

s ip stream meets the wall. In the extreme case the point of inflection seems 



46 


WALKER BLEAKNEY 


to develop into a much higher order of contact, the compression becomes a 
shock, and a second slip stream is visible, as shown in Fig. 5. 

For some time the fact has been emphasized that the three-shock theory 
failed for weak shocks rather than strong. Perhaps a more significant way 
of stating the same fact is that the theory fails for curved shocks. Some 
effects of curvature and its derivatives have been worked out by Taub [31], 
and these may throw some light on the local theory. 
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ON A NONLINEAR DIFFERENTIAL EQUATION IN HYDRAULICS 

BY 

N. W. McLACHLAN 

1. Introduction. The equation 

(!) u + 2k\u\u + au = d ( cl,k > 0,d ^ 0 ) 

occurs in connection with a hydroelectric power system. Referring to Fig. 
1, water is supplied from a reservoir via a tunnel and pipe line to turbines which 
drive electrical generators. The rate of flow of water to the turbines is con¬ 
trolled by valves at the lower end of the pipe line, which may be opened or 
closed as desired. If one or more valves are closed, in order to avoid the 
occurrence of a large force in the system due to change in momentum of the 
surplus water, a surge tank is situated at the lower end of the tunnel. On 
closing one or more valves the water level in the tank rises above its value for 
steady flow, and this restricts the force to a safe maximum. When one or more 



Fig. 1. Schematic diagram for hydroelectric system having surge tank. 


valves are opened, water flows from the tank, thereby augmenting the supply 
from the reservoir. In either case, ?/, the velocity of the water in the tunnel, 
ultimately attains a substantially steady value u = d/a. Any change in the 
rate of flow of water to the turbines is accompanied by a surge in the tank 
tunnel, and reservoir, which combination may be regarded as a form of enor¬ 
mous U tube having a nonuniform cross section. The surface area of the 
water in the reservoir far exceeds that in the tank, and so the level in the 
foimei is substantially constant during a surge. Owing to friction between 
the water and the internal surface of the U, loss occurs which varies as u 2 
Since the friction opposes the motion of the water, irrespective of its direction 
the damping coefficient is always positive and is represented by 2k\u\ The 
term au arises in virtue of the gravitational restoring force when the water 
level in the tank differs from its equilibrium value. For a tunnel of given 

49 




M N. W. McLACHLAN 

dimensions, a cc 1/A„ where A, (assumed constant) is the cross-sectional area 

of the surge tank, d is proportional to Q, the rate of water flow (cubic feet 

per second) to the turbines. Thus d = 0 signifies that all the valves to the 
turbines are closed. 

2. Alternative form of equation (1). Write u = z + c, where c = d/a, and 

(1) becomes 


(2) x + 2k\x + c\x + ax = 0. 

If the supply of water to the turbines is discontinued, d = 0 so c = 0, and (2) 
becomes 

(3) x + 2k\x\x + ax = 0 . 

It may be shown for initial conditions x(0) = x 0 , x(0) = x h that x and x-+0 
as t —-> + oo. 

3. Singular points of (1). Putting d = 0 and writing v = du/dt, v' u = dv/du 
we obtain 


(4) 

so 

(5) 


vv' u + 2k\u\v + au = 0, 


/ _ — (au + 2k\u\v) 

u 

V 


When u u 0, the right-hand side of (5) has the indeterminate form 0/0, so 

the origin is a singular point, and by Poincares criteria [1] it is either a center 

or a spiral point. Presence of the higher-order term 2k\u\v does not permit 

discrimination between these two. However, in virtue of the last sentence 

in Sec. 2, the origin is a stable spiral point. Consequently, if d = 0, the water 

executes a damped oscillation whatever the value of k > 0. This conclusion 

may be grasped more readily if it is realized that the damping coefficient 

2k\u\ -> 0 as w-> 0. Moreover, the motion is ultimately quasi-periodic with 
approximate period 2w/aK 

When d > 0, we use (2) in the form 



= f/ _ “ (ax + 2k\x + c\v) 
dx x v 


Since x —> 0 as t —> + oo, it follows that, at some t = t x > 0, c > \x\; so in 
t\ < t < 00 , \x + c\ may be replaced by x + c. 1 Then Poincares criteria 
show that if /rc 2 < a, the point x = v = 0, that is, u = c, v = 0, is a stable 
spiral point, while if k 2 c 2 ^ a, it is a stable node. In the first case the motion 
is a damped oscillation, whereas in the second it is damped but nonoscillatory 
(aperiodic). The general forms of the (v,y)-curves (obtained by the method 
of isoclines) for the three cases are depicted in Figs. 2, 3, and 4. The numbers 
on the thin curves indicate the slope of the (y,y)-curve at the crossing points. 

1 With certain initial conditions c > \x\ for h = 0; see Sec. 4. 
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In connection with Fig. 2 it is of interest to remark that if the initial conditions 
are |u 0 | > 0, |t> 0 | = o/2k, then over the range (|u 0 |,0), we have |u| = |y| = a/ 2k, 
and so |r| is constant, and the relationship between the water velocity u and 
time is linear. When u changes sign, the (u,f)-relationship ceases to be linear. 



4 Solution of (2). If K is such that the damping is not too large, we may 

employ the method of slowly varying amplitude and phase [2], Then in the 
first approximation, the form of the solution is 



w A(<) sin [oot -f- ^>(/)], 


A(/) being the amplitude and v {t) the “phase” angle at any time t ^ 0 that 

witi! w v 6 Cti T ° f * Wh ° Se mte 0f variation is «mall in comparison 
with w. Writing a — w 2 , a )t = \f/ } we have 





sin \p + c\ cos 2 \p d\p. 


u?fo c “u,To"o X7^ ( <{ C - d/a “ w 6 * < 2 >«- < 
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Fig. 3. Integral curves for u + 10|w|u + 100a = 100, that is, x -f 10|x + l\x + lOOx = 0, 
with u = x + 1. Oscillatory case with d = 100 and k 2 c 2 < a. u = 1, v = 0, is a stable 
spiral point. 

Case 1. Here |A sin \p + c\ may be replaced by (A sin \p + c), and since 
the integral corresponding to A sin \p in (8) is zero, we obtain 


(9) 

and so 




cos 2 \p d\f/ = — kcA, 


(10) A(t) = A 0 e-‘ cl , 

where A 0 is the amplitude at t = 0. Also we have 


(ID 

so that 





sin \p + c) cos \J/ sin \p d\f/ 



(12) <p = <po , an arbitrary phase angle. 

Writing <p 0 = (7 t/2) — </> 0 , and substituting from (10), (12) into (7), the formal 
solution of (2) is 

(13) x = A 0 e~ KCl sin (wt + vo) = A 0 e~ KCt cos (cot — fa). 
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Hence, finally, the velocity of the water in the tunnel in the first approxima¬ 
tion is 

J 

(14) u = x + c = .4 0 e - "‘ cos (ut — <p 0 ) H —> 

a 


A o and </> 0 being determinable from the initial conditions, i.e., those which 
obtain at the beginning of a surge after d has been altered. 

Initial conditions. Suppose that the steady rate of flow of water to the 
turbines is reduced from Q x to Q 2 by closing a valve . 2 Then as t —► -0, 
u = c = di/a, but for t > 0 , c = d 2 /a, and we have to determine the corre¬ 
sponding value of u. The initial conditions are w(0) = d x /a, w(O) = 0, with 
c = di/a in (14). Thus 


(15) 

so that 

(16) 


A 0 e~* dlt/a cos (c ot — <po) 




A 0 cos <t> o 



Differentiating (14) and putting t = 0 gives 


(17) tan </>o = ^ and cos 0 O = ^ j 

Substituting for cos 0 O into ( 1 G) yields 



Hence from (14) and (18) the velocity of the water in the tunnel during the 
surge is 





COS (c ot — 0o) -J— 



where <t> 0 = tan~ 1 (sd 2 /wa) . From (19) it is seen that as t -> + °o, u -> d 2 /a, the 
steady velocity after subsidence of the surge. By hypothesis c = d 2 /a ^ A 0 - 
so we must have di/{ 1 + [1 + ( K d 2 / a><z) 2 ]*} < d 2 < d x . 

Suppose that the steady rate of water flow to the turbines is increased from 
Q i t0 ^ 2 - Then < d *> and if d 2 /a &; A 0 , (19) is valid and we have 







e Kdlt/a cos (a )t — 0 0 ) + 



As before, when t — + =c, u - d 2 /a, the final unidirectional velocity of the 
water in the tunnel. 


2 lhe time taken for the change from Qi to Q 2 is assumed small 


compared with 2 tt/co [3]. 
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Fig. 4. Integral curves for u + 10|u|u + 100a = 200, that is, x + 10|x + 2|x + lOOx = 0. 
Nonoscillatory case with d = 200 and k s c 2 = a. u = 2, v = 0 is a stable node. 


6 . Case 2. |A| > c = d/a, 0 ^ t < ti. In this instance evaluation of the 

integral in ( 8 ) yields 


( 21 ) 


A = — 


kA f/ 2 (A 2 - c 2 ) 1(2.4 2 + c 2 ) 


n( 


3A 


+ 2c sin 




and as in Sec. 4 we assume that the damping is not too large, 
as a first approximation ( 21 ) gives 


If c/\A\ ^ 0.4, 


( 22 ) 
so that 

(23) 


A = 


- (£) 


4 tct 

A~ l = 5 -h B, a const. 

07T 


If A 
(24) 


= Ao when t = 0, B = 1/A 0 and 




1 + ( 4:A 0 Kt/37T ) 


We find also that p = <p 0 , an arbitrary phase angle. Hence by (24) and (7), 
with <p Q = (tt/2) — <f> o, we get 
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(25) 


U = X + C = 


An 


1 + 7 * 


cos (a )t — <po) + c , 


where 7 = 4 Aok/ 37 t, and provided that A 0 /( 1 + yt) > c = d/a , 0 ^ t < t x . 

When the water supply to the turbines is reduced to a comparatively small 
value so that di » d 2 , the initial conditions are ^( 0 ) = d x [a , zi( 0 ) = 0 , and 
we take c = d 2 /a. Thus, 

(d 1 — d 2 )/a 


(26) 

and 

(27) 


A 0 — 


{1 - [Mdi - d 2 )/3™a] 2 fij 


u = 


0 


1 + yt 


cos ( o)t — tan -1 - 


a + 

to / a 


provided that 4 0 /(l + tO > da/a, and 0 ^ t < t l = ( 37 r/ 4 K)(a/d 2 - 1/4 0 ) 
When £ > <i, the solution has the form at (14). viz., 


(28) 


U = Aie~ KCt COS (o;£ — </)i) + C, 


where Ai and </>i may be determined from the values of u and u at t 
Then by (25) 


— 1 1 


(29) 
and 

(30) 


u(t x ) = 


An 


1+7*1 


COS Mi - 0 O ) + c, 


MO = - 


A o 


1 + 7*i L 


7 COS Mi - </>o) . . / J 

- x ^ -b co Sin Mi - <t > 0 ) 


From (28), (29), and (30), we find that 


(31) 

and 

(32) 


^ _ Aoe* rtl cos Mi — 0 O ) 

1 1 + 7*1 COS Mi — 0!)' 


<t> i = co^i + tan -1 


KC 


Leo co(l + 7 ti) 


— tan Mi — </>„) 


I- 4 1 > c = 0. < = 0. Suppose that when the velocity is di/a, the water 
supply to the turbines is stopped completely. The initial conditions are now 

■u( 0 ) = di/a, u( 0 ) = 0 , with c = 0 , and the required results may be obtained 
from (26), (27) by putting d 2 = 0. Thus 


(33) 


°“ a/ 


1 - 


4kc/i 
37 tc oa 


and the velocity of the water in the tunnel is 


(34) 


u = 


A o 


1+7* 


cos 


^cot — 


tan 


—i 



Comparison of (34) with (19) shows that the decay of the oscillation is more 
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rapid in the latter case, owing to the respective damping factors being 1/(1 -f 
yt) and tr*#*. A flow of water entails d 2 > 0, so that in (2) the respective 

damping coefficients for tf 2 > 0 and d 2 = 0 are 2k\x + c| and 2k\x\, the former 
being greater than the latter. 

6. The equation y + 2 K yy + ay = 0 (a,* > 0). This is (1) with d = 0 

and y written for \y\, so that the damping coefficient changes sign with y. 

If the equation is that for a dynamical or other physical system, when y > 0 

some of the energy in the system is dissipated, whereas when y < 0 energy is 

supplied to the system from an external source. In the case of an electrical 

circuit, there would be a negative-resistance effect as with an electron-tube 
oscillator. 

Taking v = dy/dt , v' = dv/dy , the above equation may be written 
( 35 ) vdv + 2 Ky(b + v) dy = 0 

with b = a /2k. Integrating (35) gives 

v ~ & In (6 + v) + kij 2 = A, a const., 
provided v > —b. When v < — b, the appropriate result is 

v b l n (b + y)] + Ky 2 = B, a const. 

If the initial conditions for (36) are y = y 0 > 0, v = v 0 = 0, we obtain 

(38) v — b ln — - + K(y 2 — y 2 ) = 0, 
and 

(39) V = ± (-v + b ln ) + ygj . 

so the (y, 2 /)-curves are symmetrical about the y-axis, when v > — b. Also 

the points y = ±?/o, v = 0 are extrema on the y-axis, while on the y-axis, 

where y = 0, there are two values, one positive and the other negative (> —5), 
which satisfy 

(40) v-b In - K yl = 0. 

I hese values, of which the positive one is the greater numerically, are extrema 
on the y-axis. 

From (35) 

(41) v' = — y(a + 2kv)/v, 

and the origin y = v = 0 is either a center or a spiral point. Equation (41) 
shows that the (v,y )~curves intersect the axes orthogonally. If for any v the 
sign of y is changed, so also is that of v'. Also, the slopes of the (v,y)-curves 
are negative in the first and third quadrants, positive in the second and fourth 
quadrants, while by (39) the curves are symmetrical about the ti-axis. From 
these deductions together with the existence of extrema on the two axes, we 
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conclude that the origin is a center, which entails stable periodic motion, 
provided v > — 6 . 

Formula (38) may be written 
(42) v — b In ^1 + ^ + Ky 2 = Ky^ 

If \v/b\ < 1 is small enough, (42) takes the approximate form 



which is the equation for an ellipse with the origin as center and semiaxes y 0) ab/ 0 . 

The (t>,?/)-curves for (37) entail v and t*o < — b. Thus with initial conditions 
V = 2 /o, v = t’o < — 6 , we have 

(44) v = v 0 -b In + x(y* - y*), 

so that 





(v 0 — v) — b 111 V ° 

0 + V 



which indicates that when v < — b, the (y,y)-curves are symmetrical about the 
n-axis. Also, it follows from (44) that, with ti 0 < — b, v— > — °o as y — > —<», 

whatever the value of y 0 . Hence, under the foregoing initial conditions, the 
motion of a physical system which obeys the differential equation 


V + 2 Kiyy + ay = 0 (a,* > 0 ) 

is unstable, and y — » — <w as t—> + °c. 

The forms of the (v,y)-curves corresponding to the stable (v > -b) and 
unstable (v < —b) cases are illustrated in Fig. 5. 

7. Physical considerations. In the schematic diagram, Fig. 6, s is a uniform 
coil spring of negligible mass, whose force-displacement relationship is / = ay, 
a > 0 . m is a unit mass constrained by frictionless guides to move along the 
axis of the spring. The damping device is such that when m is moving to the 
right 3 the damping coefficient is 2 K y > 0 , while for motion to the left it is 
2ki/ < 0. Let m be displaced to the right a distance yo, thereby extending s 
by this amount. On m being released, its velocity will increase until the 
central position y = 0 is reached. During the corresponding time interval 

ro, there will be energy loss 2k yv 2 dt, while, during the succeeding interval 

ro, this energy will have been restored; so m reaches the extreme position 
y = -y 0 . On the journey from y = -y 0 to 0, 2 xy < 0, so that energy will 
be supplied to the system over a time interval n, say. Consequently the 
velocity of m in the central position will exceed that when it has moved there 
from y = y 0 . The kinetic energy of m is now greater than the potential 


3 Of the central position. 
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Fig. 5. Integral curves (y,y) for y -j- 10 yy 100 y = 0, showing stable (above the sepa- 

ratrix) and unstable regions (below the separatrix). The origin is a center. —a/2/c = — 10 . 



Fig. 6 . Schematic diagram for mass-spring system having positive-negative damping 
device. 
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energy at either y — ±y 0 by an amount 2k ' yv 2 dt. During the succeeding 

interval n, the damping device causes dissipation, and the whole of the energy 
gained in the previous n interval has been lost when m reaches y — y 0 . 

If y o is large enough, the (i>,y)-curve in the third and fourth quadrants will 
be sensibly rectangular, the velocity over most of the path being approximately 

v ~ — a I 2 k. lhus 2r 0 ~4Ki/o/a. which is the approximate time of travel 
from y 0 to —y 0 , so 2r 0 —> * as (/<>—> x . 



?“• 7 ‘ ®°! Uti0n Cu ™ s (»,<), for y + 10 yy + 100,/ = 0. Solid curve is for initial con 
ditions i/(0) — Vd — 3, i/(0) = 0; for broken curve y( 0) = y a = 0.1, y(0) = 0. 


In the first and second quadrants, so long as v is not too small the loaa 
rithmic term may be neglected in (38), and we get 


(40) 


K (yn - y-) = 


= xy'o i - ( 


y_ 

2 /o 


The time interval corresponding to the movement of m over any range (- y l>yi ) 

where l^l < |j/ 0 |, for which (46) is valid, is given by t = 2 f v ‘ dy/v, and t 

decreases with increase in y 0 . The time intervals for the ranges (~v 0 - y ) 
and (y h y 0 ) are correspondingly small. Hence the periodic time 2 (r„’+r\ 
increases indefinitely with increase in y 0 . For large finite values of the latter 

form of a “"■ tooih wave i,ith roi '" d ' d » 

To explain the (v,y)- curves below the line e = -a/2 K , consider the initial 
conditions to be such that when y = y 0> m has a velocity „ 0 < - a / 2 ” in 

the range y < 0, the control ay is inadequate to prevent v increasing (negatively) 
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under the influence of negative damping, so that both y and v —* — oo as 

t — >-)-oo. 

8. Approximate solution of y + 2k yij + ay = 0,y small, •; > -a/ 2k. If K is 

relatively small, the method of perturbation may be used, since the solution 
is periodic. Thus we assume that 


( 47 ) y = y 0 + Kyi + K 2 y 2 + • • • , 

and 


(48) a = Wq -f- kwj + + • • • , 

where y 0 , y i, y 2 , • • • are twice differentiable functions to be determined, while 
«o, «!,••• are evaluable constants, in virtue of the periodic nature of the 
solution. Proceeding in known manner and writing \p 0 = co 0 <, to the second 
order in k, we find that, with initial conditions 7/(0) = y„, 7/(0) = 0, 

(49) y = 7/o ^1 - cos ^ (2 sin - sin 2<Po) 

2 3 

+ (32 cos 2^o — 9 cos 3\po), 


where c 05 = a — (irt/o/3), and we assume that k 2 ^/3 <<C a. 

For initial conditions y( 0) = 0, ?/(0) = <a Q y 0 , to order 2 in k the solution is 





sin \f /0 + 


“Vo 

3cOo 




sin 2\f/o 



k 2 Vo 

8 


sin 3^o, 


with col as above. 

In these instances the (v,y )-curves would be approximately elliptical in 
shape with the origin as center, particularly so if a/2* were large. The periodic 
time of the motion is 



2tt ^ 27 t (1 + K 2 yl/§a) 
O)o oft 


so r increases with increase in y 0 . This feature persists however large y 0 may 
be, as was demonstrated in Sec. 7. 

9. Consideration of y + 2 Kyy ft- ay = d > 0. Writing y = x + c, we 
obtain 


(52) x + 2 k(x + c)x + ax = 0, 

provided c = d/a. The paragraph below (6) shows that if k 2 c 2 < a, the point 
y = c, v = 0 is a stable spiral point, while if k 2 c 2 ^ a, it is a stable node which 
the (v,y)-c\irve enters as t — > + °o. Thus the respective motions corresponding 
to these conditions are a damped oscillation and aperiodicity (absence of 
oscillation). When k 2 c 2 « a —the values of ?/(0), i/(0), and k not being too 
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large — in the first approximation the solution is given by (14), so that with y 
for u we have 

(53) y = A 0 e~ KCt cos (cot — <p 0 ) + -• 

a 

If the initial conditions are y = y 0f v = 0, then 

(54) y = (y°- f) [i + (^) ] e-'“ cos (U - 4 , 0 ) + i 
with 0o = tan" 1 ( K d/ua). 
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ACOUSTIC RADIATION PRESSURE ON A CIRCULAR DISK* 


BY 


HAROLD LEVINE 

1. Introduction. The radiation pressure experienced by material objects 

in a wave field represents an important acoustical effect. It is therefore of 

interest to have the related theory available for a variety of experimental 
conditions. 

A calculation of radiation pressure on a plane circular disk of rigid material 
by King [1] furnishes the inspiration for this paper. The disk is exposed to 
small-amplitude plane progressive or stationary harmonic waves with equi- 
phase surfaces parallel to its plane and performs oscillations along the normal 
direction. In the analysis of King, both inertial and diffraction effects are 
taken into account approximately, the latter only for wavelengths large com¬ 
pared with the disk radius. We present here another method of calculating 
the radiation pressure which is capable of high accuracy for all wavelengths. 

King’s procedure employs integral representations of the wave function 
and leads to a pair of dual integral equations, one in the domain of the disk 
and the other in the remainder of its plane. An approximate solution of these 
equations is effected by appealing to the known solution of related frequency- 
independent or static integral equations. This process furnishes the first few 
terms in the low-frequency expansion of the solution, with succeeding terms 
obtainable after tedious computations. Estimates of the radiation pressure 
are thereby restricted to low frequencies or long wavelengths. 

The scheme proposed here involves a single integral equation which deter¬ 
mines the difference in velocity potential on opposite disk faces. From the 
integral equation and its connection with the normal radiation pressure, we 
construct an expression which is stationary for variations about the correct 
solution and has a stationary value equal to that of the radiation pressure. 
The stationary form has the practical merit of subordinating difficulties asso¬ 
ciated with the construction of an exact or generally accurate solution to the 

integral equation. We shall illustrate this approach by extending the results 
ot King to include the high-frequency limit. 

2. Formulation. At equilibrium we locate the disk in the plane z = 0. 
en the incident radiation field has axial symmetry with respect to z, a 

circular disk will carry out oscillations perpendicular to its plane, with no 
tendency for angular deflections. If *(r) denote the velocity potential of a 
held possessing harmonic time dependence exp (-iwt), the conditions of the 
problem are met by constructing a solution of the wave equation 


versTtv S J 0 ^l f VaS • Per !° rm ! d °“ a Nati ° nal Bureau of Standards contract with the Uni- 
Research, U.s'Kavy ° S gelpS ’ and Was s P on sored (in part) by the Office of Naval 
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W (V 2 + k*)*(T) =0, k = - 

c 

(c denotes the sound velocity) which incorporates the prescribed form of inci- 
dent radiation and has normal derivative 



d , / N 


on the disk, f being the disk velocity. The latter quantity is taken from the 
dynamical equation for the disk, whose motion is ascribed to the first-order 
pressure difference on opposite faces. Employing the scalar Green’s function 



G(r,r') 


exp (ik\i — r'|) 
47r|r — r'l 


G(r',r) 


) 


which satisfies the inhomogeneous differential equation 
W (V 2 + A; 2 )G(r,r') = — 5(r - r'), 


a representation of the velocity potential deduced from Green’s integral 
theorem is 


(5) ^(r) = ^"‘(r) + J[G(r',r) d n ^{r') - *(r') d n tf( r ',r)] ds', 

wheie the integral extends over both disk faces. Since the derivatives in 

(5) lefer to the inward normal and the disk moves rigidly, d n >\// has equal and 

opposite values for any pair of points on the two faces. Then the first term 
of the integral vanishes, giving 

( 6 ) Hr) = H nc (r) - J a \k(r') d S 'G(x',y',0,T ) els', 
where 


( ? ) = H*,y,z -> 0-) - i{x,y,z -> 0+) (x,y in A) 

and A denotes the disk area. 

With the help of the Green's function Fourier-integral representation 


(8) 

f’° exp {i[k x (x - x') + k u (y - y') + \/k 2 - k\ - k* \z - z'\]} (jj . 

J “ y/k 2 — kl — kl 

and the Dirac delta-function relations 



5(x - x') b{y - y') = (2x) 2 j _ exp j i[k x (x — x') + k y {y 

I\f(x',y') K* - *') S(y - l /') dx' dy' = f(x,y), 


it follows from (6) that 



i{x,y,z -> 0-) = \p inc (x,y,z -> 0-) + |¥(r), 
Hx,y,z -» 0+) = i ine (x,y,z -* 0+) - |4-(r). 


y')] \ dkzdk,, 


(x,y in A), 
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On subtracting the latter equations and noting that 

t ine (x,y,z -* 0 -) = \l/ ine (x,y,z -> 0 + ) = t ine (x,y, 0 ), 
we obtain the identity 

*(r) = *(r), 

whereas addition yields the new result 

( u ) W*#,* -> 0 -) + $(x,y,z -> 0 + ) = 2\l/ ine (x,y,0). 

For the progressive incident wave 

(12) \l/ inc (r) = e ikz , 
we satisfy the boundary condition (2) if 

( 13 ) -t = ik - f *(r')K(r',r) ds' (r on A) 


04) 


K( r,r') = d,d z ’G(x,yfl,x',y'fl) = K(r',r). 


lo write a dynamical equation for the disk requires information about the 
pressure in the medium. It is well known that with inclusion of quadratic 
or second-order terms in the velocity potential the deviation of pressure p 
from its equilibrium value p 0 takes the form [1] 



where p 0 is the equilibrium density of the medium. Using only the linear term 
in (15), the equation of motion for a disk of mass m 0 admits the first integral 



Ihus, by eliminating f between (13) and (16), we arrive at an integral equation 
for the function 'k(r), 




*(r')K(f,r) ds' 


Po_ 

m 0 



Sl'(r) ds = ik y 


(r on A)y 


which is specialized to the rigidly fixed disk if m 0 -> «. 

We next inquire after the time-average radiation pressure on the disk in the 
noima or z direction, the linear term in (15) evidently does not contribute 
to the result, and likewise for the final term, since 

I a U^(w-> o-))* - 0 + ))*]ds 


in 

be 


= 2 f A diC‘(r) -f ^(r) d„C'(r)] ds = 0, 

consequence of (2), (7), ( 11 ), aIld (12) . Finally, the normal pressure 
comes 
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0-))2 - (4,(x,y,z-* 0+)) 2 ]ds 


or, on taking the time average, 


< i9 > p - 4’X m *• 

A different characterization of the time-average radiation pressure is 
obtained by employing the integral equation (17) for *(r). On multiplication 
of the latter equation by ¥(r) and subsequent integration over A, there results 

(20) L *W^ r ’ r ')*( r ') ds ds' - ^ ¥(r) dsj = ilc *(r) ds, 

whence 

(21) II ^ (r)Z(r ’ r ')^( r ') ds ds ’ - ^ ( J A *(r) dsj - 2ik *(r) ds 

— —ik J 'k(r) ds. 

One easily verifies that the left-hand member of (21) is stationary for variations 

of *(r) about the solution to the integral equation (17), with the right-hand 

member as its stationary value. Moreover, by utilizing the relation (19), 

we deduce that the time-average radiation pressure is the stationary value of 
the functional 


(22) P(=P(*)) = -ipofclm ^ 


'k(r)if(r,r , )'k(r') ds ds' 



this feature imparts practical significance to (22) as a means of approximating 

the radiation pressure. In the latter connection, it is useful to have another 
form of (22), 


(23) P = - 1 Po A) 3 Im 



'P'(r) ds 



ds ds' - ^ 
a m 0 



^(r) ds 



which is independent of the scale for \k(r). 

Before applying (23) to the specific case of a circular disk, we may remark 
on some other general features. Thus, when the disk has infinite inertia 
(m 0 = 00 ), the radiation pressure and plane-wave scattering cross section at 
normal incidence, <r ac , of the disk are simply related by 
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P = vpok 2 <r 


eey 



or 

(25) 

where 





is the time-average energy density for the incident wave (12). The relation 
(24) is substantiated on exhibiting a stationary scale independent form of the 
cross section in terms of 4>(r), for details of which we refer to another paper [2]. 

In the limit of extremely short wavelengths, or very high frequencies, the 
normal scattering cross section for a rigid disk approaches a value of twice its 
area A , and consequently 

(27) P ~ 2W inc A, (fc— oo). 

This result can be understood in terms of the momentum transfer which 
accompanies reflection of the incident plane wave at the disk, with the factor 
2 specific for unit reflection coefficient at a rigid surface. Even if the disk has 
a finite mass mo, the asymptotic relation (27) holds, for the first term in the 
denominator of (23) then dominates the second. 

At low frequencies, examination of the integral equation (17) shows that 
'l'(r) is wholly imaginary and (except for a scale factor) independent of fre¬ 
quency. With this knowledge, and use of the expansion 


(28) 


= *.( r,r') + 


ifc 8 

127T 


+ 


K,( r,r') = d g d a 


1 


r — r' 


z = z' *=0 


it follows by (23) that 

(29) 

where 


P 2 pok (T.c ^ 


(30) 


and 

(31) 


m\ = p o 


Ua * (r) ds ) 


j A *(r)K.(r,r')*(r') ds ds’’ 


_ k 4 ml 


(Tbc = 


127T P P 








I he inertial effect is thus isolated in a single factor and characterized by the 
quantity in i, the so-called hydrodynamic mass of the disk. 

A final remark in connection with the actual use of (23) is appropriate; 
as K(t,t') is a highly singular function, nonvanishing estimates of the radia- 
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tion pressure are obtained only for the restricted class of trial functions *(r) 
that vanish on the rim of the disk. 

3. Application In regard to a circular disk with radius R, the function 

(r) may be adequately described by linear combinations of the primitive 

orm [1 - (r 2 /E 2 )]"+J, n = 0, 1, • • • , where r denotes the radial coordinate 
re erred to the disk center. Thus, as a first approximation, we take 


(32) 


*i(r) = 


- - f )'■ 


in (23), and avail ourselves of the results in [2] to deduce that 


(33) 

with 


P l = 1 pJkR) 3 _ Iu(kR) _ 

9 [R n (kR) + i(poR 3 /m„)] 2 + [/n(/c/2)] 2 ’ 


In(a) = £ - 


and 


2x *xo + 8a* + Xa ~ (’ + J^i) [ *«)<*. 


W " 1 0 + i) / 0 ! " * - £1 - S 


J i(2 a). 


Heie J 0 , J 1 and $ 0 , aSi denote Bessel and Struve functions of order zero and 1, 

respectively. Using well-known properties of these functions, it follows that, 
for small values of kR , 


whence 


1 u(UR) = ~ {hRY, R u (kR) = g, 


(34) 


Pi 


27tt 


P0 (kR) 


[1 + %pa(R 3 /m 0 )] 


(kR « 1) 


A comparison of (29) and (34) establishes 


(35) 


mi = |p 0 /? 3 , 


as the hydrodynamic mass for the circular disk with radius R. 

The result (34) also emerges from King’s calculations, which, however, 
disagree slightly in the succeeding frequency correction. Thus, with two terms 
for the estimate of Rn(kR), kR <3C 1, equation (23) yields 


(36) 


p, = A 0n( k R) 6 1 +^ kR y- 

27tt M ; [1 + %p 0 (R 3 /m 0 )(l + i(W?) 2 )] 2 ’ 


while King gives the (incorrect) numerical coefficient 4 instead of -£%. 

At very high frequencies, where King’s procedure becomes impractical, the 
asymptotic behavior of Pi, 


(37) 


Pi W»(kR ) 2 , 


(kR » 1), 
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is easily derived, since then 


In(kR) 



Rn(kR) 0. 


To ascertain the approximate nature of (37), we rewrite it in the form 

(38) Pi^frAK&R 2 ), 

which permits ready comparison with (26), (27). This shows that (38) devi¬ 
ates from the correct result by a numerical factor in the scattering cross section, 
having ^ in place of 2. For intermediate values of kR, numerical values of 
Pi are obtained from (23) without difficulty. 

An improvement in these results is achieved by use of the trial function 



where c is a variational parameter, that can be specified by the stationary 
condition 


dP 2 

dc 



In this way we derive an estimate P 2 of the radiation pressure which agrees 
precisely at low frequencies with Pi and is generally more accurate at higher 
frequencies, attaining asymptotically the form 

( 4 °) P 2 ^ W^iVbrR 2 ) = UirpoikR)*, (kR » 1). 

I he case of excitation by stationary plane waves is treated analogously. 
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INFINITE MATRICES ASSOCIATED WITH A 

DIFFRACTION PROBLEM 1 

BY 


W. MAGNUS 


The problem of diffraction of a plane scalar (acoustical) wave by a circular 
aperture in an infinite plane screen leads to the following boundary-value 
problem: Find a solution u of V 2 w + k 2 u = 0 (with a constant k) with the 
following properties: If z,p,d are cylindrical coordinates, u satisfies the differ¬ 
ential equation everywhere except for z = 0 and p ^ a > 0. For 2 = 0, 
P > a } u vanishes, and for z = 0, 0 ^ p < a, du/dz assumes a given constant 
value Vo ^ 0. At infinity, u satisfies a Sommerfeld-radiation condition. Then 
u = u(p,z,0) is uniquely determined and represents the diffracted field. 

Because of the symmetry of the boundary conditions, u does not depend on 
0. In the diffracting circular aperture, z'.e., for 2 = 0, 0 ^ p ^ a, u = $(p), 
where $ depends on p only. It has been shown by Levine and Schwinger [1] 
that u can be determined completely if 3>(p) or even if C 0 ^(p) is known where 
Co denotes an undetermined constant factor different from zero. 

Expanding <£(p) in a series 


ao 



m = 0 


and applying their “ variational method,” Levine and Schwinger show that the 
Xm s satisfy a system of inhomogeneous linear equations 


^ ln,m(P)x m = 7? n , (tt = 0,1,2, * * •), 

m = 0 

where the l n ,m’s are explicitly known power series in @ = ka /2 and where 
Vn = l/(zi + f). Since the quantity 

( 3 ) t* = y — x 7 -3 

m % m + * 

determines the “transmission coefficient” of the diffracted wave which is used 
in deriving (2), only those solutions of (2) will be considered for which the 
series in (3) converges. Then the series in (1) converges for 0 < p < a; since 
3>(p) must be finite and continuous for p = 0, we also assume that 


1 Abstract. This work was performed at the Washington Square College of Arts and 
Science, New York University, and was supported in part by Contract AF-19(122)-42 
with the U.S. Air Force through the sponsorship of the Geophysics Research Division Air 
Force Cambridge Research Center, Air Materiel Command. 
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exists. 


The l n . m ’s are given by 



00 

^ ln,m0 q > 

5 = 0 

where for q = 2 p, p = 0 , 1 , 2 , • • • , 


^( 2 p) _ ( l) p 7 r^ r(m -f~ %)T(n -j~ -§-)r (m -f- n -f~ 2 p -|- 1 ) 

4p! T(w + p + l)r(n + p + l)r(m + n + p + 


and for q = 2p + 3, 


j( 2 p+ 3 ) _ l) p 7r^ r(m -f- |-)r(n ~t~ |-)r(m -f- n -f- 2p -f~ 4) 

4F(p + f) T(m + p + f)r(n + p + f)r(ra + n + p + 4) 

For p = 1, /nVm = 0. Denoting by L (p) the matrix with the element V£ m in its 
(n + l)st row and (m + l)st column (n,ra = 0,1,2, • • •), it can be shown that 
for p = 0, 1, 2, 3, • • • 

L (p) = L ( 0 ) aS ( p) , 

where $ (0) is the identity and where the matrices S (p) are known explicitly. 
For odd values of p, all the elements in £ (p) except those in the first (p - l)/2 
rows are zero, and for even values of p, the element in the (1 + n)th row and 
(1 + m)th column of *S (p) is zero if 2n < p + 2m. The £ (p), s are bounded 
matrices, and estimates for their upper bounds can be found. From this it 
follows: 

If the linear equations 

oo 

( 5 ) 2 l ^ 0) = Vn 

m = 0 

have a solution oi 0) satisfying the conditions implied by (3) and (4), then the 
aim's can be represented by power series in 0 , 



which converge at least in a certain range 0 :g ^ /3 0 . The x (p) ’s can be 
derived from the x^ 0), s by finite-recurrence relations. In the particular case 
where y] m = l/(m + f) we have a:£ 0) = 8 / 7 r, x^ = 0 for m > 0 , and at most 
the first p + 1 of the numbers x {p) are different from zero. They can be deter¬ 
mined by using only the first p + 1 of equations ( 2 ) and neglecting all terms 
involving x p+h Xp+ 2 , * * * and all terms involving 0 p+1 , fi p+2 , • • • . The 
aim's satisfy the conditions implied by (3) and (4). 
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The factorization of the L (p) ’s into a product L (0) <S (p> shows that the solution 

of (2) depends largely on the inversion of the matrix L (0) , which represents 

the static case k = 0 (or /3 = 0). For the limiting case /3 —» Levine and 

Schwinger have shown that the variational method leads to a system of linear 
equations 


00 


( 6 ) 


2 = c 2Vn , 


m — 0 


where C 2 is a constant and where 


(7) 


ITJ = (» + m + 2)" 1 . 


It can be shown that both (6) and (5) are equivalent to the solution of a prob¬ 
lem of moments. Consequently, it may happen that (5) or (6) does not have 
any solutions and in particular does not have a solution satisfying the conditions 
implied by (3) and (4). But even if (5) or (G) has a solution for a particular 
choice of the s, no solution will exist if only a finite number of these yjs is 
changed. As a matter of fact, (6) does not have any solution at all if 


Vn = 


c 


3 


n + 7>- 


Nevertheless, the quantity T* can still be determined and even computed from 
(5), and since it is the aim of the variational method to calculate the trans¬ 
mission coefficient (that is, T*), the method still works in this case. To prove 

this, (5) and (7) arc shown to arise from certain integral equations. 

If we define 


ao 


then the integral equation 

( 8 ) 


f(v) = £ h(z) = c 2 2 

171 = 0 n = 0 


J Q l - vz)~ l dv = h(z) 


has a solution /(») which is analytic for 0 ^ | v | < 1 if ( 6 ) has a solution for 
which 1 and C i in (3) and (4) exist. But (8) may have a solution which is 
not analytic for |»| < 1 and which therefore does not lead to a solution of (6). 

is is due to the fact that the integral operator in (8) is an extension of the 
matrix operator in (6). 

In particular, let 


Vn = 


n +4 


or 


/. 


—i 


dv. 


casc, C 2 V and although there do not exist any x (co) ’s satisfying 

(6), the quantity T* can still be defined by writing satisfying 
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f(v)v * dv = T*. 


This gives T* — C 2 . Actually T* can be computed approximately by solving 

a finite number of the equations in (6) for an equal number of unknowns and 

by introducing the “approximate” values of the x^’s thus obtained into the 
left-hand side of (9). 
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ON THE COUPLING OF TWO HALF PLANES 1 


BY 


ALBERT E. HEINS AND HERMAN FESHBACH 

1. Introduction. This paper is a continuation of some earlier work which 
dealt with the coupling of two acoustical ducts of rectangular cross section 
whose walls were composed of different acoustical materials [1]. We shall be 
concerned here with the effect of a plane wave incident upon two semi-infinite 
half planes joined along a line x = 0. Each of these half planes is assumed to 
have an acoustical property which may be described by a complex parameter, 
the admittance of the material. An approximate solution to this problem 
has been given by Morse and Bolt [2] for the special case of normal incidence; 



however, these workers did not find the scattered field, nor did they indicate 
the nature of the acoustic potential along the two half planes. We shall solve 
this problem by an appeal to the theory of integral equations. This formula¬ 
tion has been indicated in [2], where Y\ and Y 2 , the admittance parameters, 
are considered functions of a coordinate. We shall consider the admittance 
parameters as constants, and hence the integral equation is one of the Wiener- 
Ilopf type [3]. The complex Fourier transform is the appropriate tool with 
which to solve this problem, and we shall see that our methods lean heavily 
on the theory of functions of a complex variable. 

2. The Formulation of the Problem. Let us consider the half space y < 0 

whose boundary y = 0 is made of two different acoustical materials. A plane 

wave whose propagation normal makes an angle 6' with respect to the positive 

x-axis is incident upon this structure. (See Fig. 1.) The steady-state velocity 

potential 0(x,?/), with the time dependence suppressed, satisfies the two- 
dimensional wave equation [2]. 



<t>zx + <t> vy + k 2 <f> = 0. 


1 This paper is based on research 
113 between the Office of Ordnance 
Institute of Technology. 


conducted in part under Contract No. DA-36-061-Ord- 
Research of the Department of the Army and Carnegie 
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^ e shall assume for the present that the half space is slightly absorbing and 
that /v therefore possesses a small positive imaginary part. 

The acoustical surfaces y = 0, x ^ 0 may be described in terms of a complex 
parameter Y, the admittance of the surface. We shall denote the admittance 
of x < 0, y = 0 by Y i and that of x > 0, y = 0 by Y 2 . Y may be written in 
terms of its real and imaginary parts as Y = 7 — ia. 7 is the conductance 
and a is the susceptance of the material. 7 is nonnegative, and if it does not 
vanish, absorption results. For such acoustical surfaces as we have consid¬ 
ered, we shall employ the boundary conditions [ 2 ], 

4>y = iwpYi(l>, (y = 0 , x < 0 ), 

and 

0i/ = io>pY 2 <t>, (y = 0, x > 0) 

Here p is the density of the medium y < 0, and o> is the angular frequency of 
the time variation. 

It is well known that, for a half space y < 0, the scattered part of <l>(x,y), 
which we denote here by </> 8 (x,y), has the following representation: 

(4) v(x,y) = f ”" [G{x,y,x'fl)<t>' v ,(xfl) - fi)G y \x,y,x'fi)] dx', 

where G is a Green’s function associated with equation (1). Since there are 
no contributions from infinity, this has been derived subject to the conditions 
that <t > 8 as well as G satisfy the re-radiation condition of Sommerfeld for 
r = (x 2 + ?/ 2 )* —» 00. Now <f> = + </>S where </>* is the incident plane-wave 

field so that equation (4) becomes 

(5) 4>‘(x,y) = fp mx,y,x'M<t>Ax',0) - ik^x', 0)] 

- [<l>(x',0) - 0 < (* / ,O)]G lf '(a:,y,* / ,O)} dx'. 

In particular if we select a Green’s function which satisfies the boundary 
condition G y = iwpYiG at y = 0, equation (5) reduces to 

<l> 8 (z,y) = j* (fa - PJGixMyX'Mtix'fl) dx' + Ai 0 r (.r,y). 

Here </>’ = exp [i(xk x + yk u )], <f> r = exp [i{xk x — yk y )],pi = iupY i, 0 2 = iupY 2 > 
and Ai = (k v — wpYi)/(k y + copFi). Observe that <t> r (x,y) is the geomet¬ 
rically reflected plane wave. In particular, the total velocity potential has the 
representation 

<t>(x,y) = (02 - P i)G (x , y,x', 0 ) <f>(x',0) dx' + 4>*(x,y) + A^fay). 

For y = 0, we have the integral equation 

( 6 ) <t>(x, 0) = J 0 " (P 2 — 0 i)(j(.t,O,x',O) 0 (x',O) dx' + </>*(x,0) + Ai 0 r (x,O). 
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3. The Green’s Function. We shall now give an explicit form for the 
Green’s function which we employed in Sec. 2 . Let us recall that G(x,y,x',y') 
satisfies the nonhomogeneous wave equation 

(7) G xx + G uy + kH} = — <5(x’ — x')8(y — y'), 

where 5(x — x') is the Dirac delta function. More precisely, G satisfies the 
homogeneous wave equation save at the point x = x' and y = y'. At that 
point, G possesses a logarithmic singularity and hence accounts for the repre¬ 
sentation in equation (4). In addition to satisfying equation (7) we require 
two further conditions which we have employed in Sec. 2 . First we require 
that 6 satisfy the boundary condition G tJ = ioipYiG, and second we require 
that for r = (x 2 + y 2 )* —> cc we have that G is asymptotic to exp ( ikr)/A in 

so fai as its r variation is concerned. 1 hat is, G acts as a two-dimensional 
source for outward going radiation. 

^ e proceed in a formal manner with this task and indicate how one may 

verify that all the conditions imposed upon G may be satisfied. A direct 

application of the complex Fourier transform to equation ( 7 ) (for the variable x) 
gives us 

( 8 ) (r(x,y,x ,y ) ~ 2 ^ J c ex P l iz (x — %') — iay)(A cos ay' + B sin ay') dz , 

- (y < y')> 

= 2tt Jc eX ^ ~~ ~ ic *y']( A cos <*y + B sin ay) dz, 

where ^ ^ ^ ^ 


A = ---, B = — - 

a + u)pY i a(a + (cpY\) 

and a = (k 2 - z 2 )i. Here we choose the branch of a which reduces to k 

for z = 0 . It remains to describe the nature of the path C. 

The representation ( 8 ) has been derived subject to the fact that G has a 

specific asymptotic form for |x| -» =c. This implies that the integrand in ( 8 ) 

is regular in the strip |Im z\ < |Im k\ provided the denominators in A and 

B do not vanish in this strip. Subject to this condition, C is therefore a path 

in this strip. In order to cast ( 8 ) into a more useful form, certain deformations 

ot the path C are in order. Before proceeding with this task, however, let us 
note that ( 8 ) may be rewritten as 


(9) G(x,y,x’,y') = ± f ex p [iz(x - x') - iay) f e - xp 

~ T J c L 2a 


) , exp (- iay ') 

2a 


h r 


wpYi exp ( — iay') 1 , 

a(a + c opYi) dZ) ^ < y ') ; 


exp [iz(x - x’) - iay’] [ {iay) + ex P ( ~ i<*V) 

L 2 a 2a 


upYi exp ( — iay)~\ 

a(a + U pY,) dZ ’ 


(V > y’)\ 
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= { [lc[(x - xV + (y - y y]i] 

+ ^H^{k[(x-xr + (y + y , m 

_ i<ap Y 1 r exp [iz(x — x') — ia(y + y')] dz 
2tt Jc a(a + wpF i) 

The H§ ] is the Hankel function of the first kind. Incidentally the first Hankel 
function in the Green’s function accounts for the two-dimensional source in 
the half space y < 0, while the other Hankel function and the integral will 
be shown to combine with the first Hankel function so that the boundary 

conditions on the line y = 0 are 
satisfied. 

Deformationjof the path C is now 
required if we are to reduce the integral 
to a usable form. In order to do this, 
we introduce radial cuts in the 2 -plane 
(see Fig. 2). By Cauchy’s theorem, 
the path C may be deformed into the 
hyperbolic path [4 ],z = k cos (\p + tr), 
— co < r < co. The angle ^ is de¬ 
fined by 

x — x' = R cos 
y + y' = R sin 

The branch of a which we desire is 
equal to —A; sin (\p + it), and hence 
there are no critical points in the 
denominator of the integrand. More precisely, the denominator can vanish 
only if 

fci sin cosh r — fc 2 cos sinh r = copyi 
k 2 sin \f/ cosh r + ki cos \// sinh r = wpo -1 

In these equations we have a measure of k 2 small. If k 2 cos < — k\ sin 
then there is no real value of r which satisfies the first equation since y is 
positive and the left side will] always be negative in the interval — tt < \p < 

0. Hence the integral along C is the negative of the integral along the hyper¬ 
bolic arc since the contributions from the closing circular are vanished for 
x > x' when hyperbola (i) is employed. A similar situation holds for x < x' 
when hyperbola (ii) is used. It follows that 

/m\ 1 [ exp [tz(x — x') — ia(y + y')) dz 

2t rijc a(a + cpY 1 ) 

1 [” exp (ikR cosh r) dr 

~ 27t J- oe copFi — k sin + it) 

1 f * [exp (ikR cosh r)](wpFi — k sin \p c osh r) dr _ 

= “ r Jo (copFO 2 - 2*«pFi sin ^ cosh r + A: 2 (cosh 2 r - cos 2 J) 
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Because of the presence of the small positive imaginary part in A;, the integral 
in (10) is absolutely convergent for R 9 * 0. It is therefore permissible to dif¬ 
ferentiate under the integral sign. Let us denote this integral by 
Then f v — iwpYif is equal to 


i 

7r 



CO 


[exp {ikR cosh r)](u)p}'' 1 — k sin \p cosh r) 2 dr 


0 (topFi) 2 — 2 kupY 1 sin \p cosh r + fc 2 (cosh 2 r — cos 2 \p) 



1 

* JO 


00 


exp {ikR cosh t)\!/ v 


d_ 

d\k 


ojpTi — k sin \k cosh r 


L (copl'i)- — 2A'ajpFi sin \p cosh t + A ,2 (cosh 2 r — cos 2 \{/) 


dr 


provided R 5 * 0. This in turn may be rewritten as 


1 

7 r 



00 


0 


exp {ikR cosh r) (/r 


ik 2 cos 2 \p 


1 


7T 

00 



00 


sinh 2 r[exp (z7c/i cosh r)] dr 


0 (cop 1 ) 2 — 2/rcop)'i sin \p cosh r + A; 2 (cosh 2 r — cos 2 i/0 


exp {ikR cosh r) 

7T JO d\P 


opYi — k sin \f/ cos r 


_(wpKi) 2 — 2ko)pYi sin \p cosh r + A: 2 (cosh 2 r — cos 2 


dr. 


Upon integrating the second integrand by parts we find that, after some sim¬ 
plification, only the first term remains. Hence the operation — iwpYi + d/dy 
applied to the Green’s function vanishes when evaluated at either y or y' = 0. 
We have thus shown that the boundary condition at y = 0 is satisfied. In a 
similar fashion we can show that the Green’s function satisfies the homogeneous 
two-dimensional wave equation save at the points x = x' y y = 7/'. Indeed 
for y < 0, y' < 0, the last two terms in the Green’s function (9) satisfy the 
homogeneous wave equation, while the first term accounts for the unit source. 

One condition remains to be discussed— i.e., the condition of outward-going 
radiation. Again for r = {x 2 + y 2 )* —> cc the first two terms of the Green’s 
function (9) obey the requirement that they be asymptotic to exp {ikr) /r* with 
the point P , {x t i y t ) located in the finite part of the plane. The third term requires 
some further comment. Let us write cosh r = 1 + v. Then the integral 
in the Green’s function transforms into 


( 11 ) 

_ }_ f _ exp {ikRv + ikR)[upY 1 — A:(l + v) sin \p]v~^ du _ 

^ Jo (2 + v) 2 [k~{l -f- v ) 2 2 kupY 1 sin \f/{l + v) + (wpF 1) 2 — k 2 cos 2 \f/] 

In this form it is clear that (11) approaches zero as \kR\ —> 00 for Im k > 0 
as a direct consequence of the Riemann-Lebesgue lemma. But more than this 
can be said, for it is possible to give the precise asymptotic form of the integral 
as |A7?| —> co, Im k > 0. There is another lemma, due to G. N. Watson [5], 
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which states that since the integrand in (11) is regular and of bounded varia¬ 
tion lor y > 0 save for the factor ir*, we have that for 1/bLR oo I m k > 0 
(11) is asymptotic to > - , 

_L_ ex P (ikR) exp (iV/4) 

a/2tt y/kR (upY i — k sin \p) 

which is clearly of the required form when the point P'{x',y') is located in the 

finite part of the plane. We have thus verified formal derivation of the Green’s 
function. 

4. The Fourier Transform of Equation (6). In order to cast equation (6) 

into a form which is amenable to the methods of Wiener and Hopf, we write 
it as follows, ’ 


(12) 0o(x) + <t> i(x) = Mx) + (0* - 00 Mx')G(x,0,x',0) dx, 


where now 


Mx) = (3 > 0 ), 

= (x < 0), 

Mx) * 4>(X, 0), (3 < 0), 

= 0’ (x > 0), 

M x ) = (1 + Ai) exp ( ik x x), (x > q), 


(x < 0). 


Clearly then 4 n(x) is of the order 



exp [ik(x' - x)]<t>o(x') dx’ 

(*' - x)i 


as x —> oc. It follows then that <j>i(x) is dominated by the exponential 
growth, exp (k^x), where ki + ik 2 = k. Furthermore, the above integral will 
converge absolutely if 0 o (x) = 0[exp (-* fa)], x -» «, and Im S < k 2 . We 
shall assume for the present that <t> 0 (x) and <£i(x) are integrable for finite x and 
verify this assumption with the solution of equation (12). With this inte- 
grability assumption, we can say something about the regions of regularity of 
the unilateral Fourier transforms of <j> i(x) and <j> 0 (x). Since the exponential 
growth of a function is the determining factor in the location of the half 
planes of regularity (assuming integrability for finite x), we have that 


4>o(z) = / exp (— izx) 0 0 (x) dx 


is regular in the lower half plane Im z < Im 5, while 


4>i(z) = / exp (—izx)<t> i(x) dx 


is regular in the upper half plane Im z > —k 2 . 
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The function G(x, 0,0,0) has the bilateral transform i/(upY\ + a ), where the 
root of the radical a has been chosen to reduce to k for 2 = 0. In view of the 
fact that this transform has only branch points as singularities, it is regular in 
the strip — Im k < Im z < Im k. Indeed it is actually regular in the cut 
plane, with such cuts as the radial ones described in Sec. 3 as possible useful 
ones. Finally <f> 2 (x) has the unilateral transform (1 + A\)/i(z — fc x ), which 
is regular in the lower half plane Im z <k 2 cos 0'. There is then a common 
strip of analyticity for these transforms that is 


— fc 2 < Im z < k 2 cos 0' < or Im <5. 

The right side of this inequality is to be interpreted as the smaller of the two 
terms. Because of the existence of this common strip, we may apply the 
Fourier transform theorem to equation (12) to obtain 


or 

(13) 

The explicit solution of equation (12) depends on our ability to decompose 
equation (13) into two parts, one of which is regular in the upper half plane 
Im 2 > — k 2y while the other is regular in the lower half plane Im 2 < Im 8 
or Im k x . Our next task is therefore the decomposition of 


♦.a + *,(«) - J±4i, + - «*■« 


1(2 - Kx) 


a + 03 p Y1 


, / n 1 . , . a + 03 p Y -i 1 + * 1 1 

$1(2) + <*>«(2) , = T- , * • 

a + oipY! l(z K x ) 



a + 03 p Y 2 
a -\- 03pY 1 


into two such functions as we have just described. 

5. The Decomposition of K(z). The required factoring of K(z) may be 
carried out with the aid of Cauchy’s second integral theorem. In order to 
accomplish this, we note that the multiplicative decomposition of K(z) is 
essentially equivalent to the additive decomposition of the logarithm of K(z). 
From the calculational point of view, however, it is better to examine the 
additive decomposition of the derivative of the logarithm of K{z). Let us 

note first that K(z) has two critical points at ±k. We shall consider here 
however, ’ 



M(z) = 1 + 



individually. After having obtained the required decomposition of each of 
these functions, we may form their ratio to find the decomposition of K(z). 
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Within the strip indicated in Fig. (3), L(z) as well as its logarithm is regular, 
he arguments of the radicals on the branch cuts (which are now drawn parallel 
to the horizontal axis of the z-plane) are indicated in the figure. Then 

L'(z) _ 1 f wpYit dt 

L{z) 2m J c [( k* - t 2 )i + wpFdOfc* - t*)(t-z) 

= _ oopYiZ 

[(k 2 - Z 2 )i + upYi](/c 2 - z i) 

The path C is a rectangular one which includes the point z and is indented above 
the branch point at z = —k and below the branch point at z = k. The vertical 



Fig. 3 


sides of this rectangle give rise to terms which vanish as they recede to infinity 

in either direction, in view of the fact that the integrand is 0(r 3 ), |/| -> ». 

Cleaily, then, the integral along the lower side of the rectangle, for z not on C 

but inside the strip, provides us with a function regular in the upper half 

plane Im z > -k 2 , while the integral along the upper side of the rectangle for 

z not on C provides us with a function regular in the lower half plane Im z < k 2 . 

Let us denote these component paths of C by C+ and C_, respectively. Useful 

forms for the integrals along C+ and CL may be obtained by appropriate 
deformations. 

We consider the integral along C+ first. Let us introduce a closed path C" 
which passes under the left branch cut, is indented below the branch cut 
t = —k, and is closed by a semicircular arc in the lower half plane. The semi¬ 
circular portion of C' + provides a vanishing contribution in the limit as its 
radius becomes infinite. Hence the horizontal portion of C' + combines with 
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C+ in the following manner when we take into account that we have a simple pole 
as well as a branch point at t = —k: 


L' + (z) 

L+(z) 


2tti j c+ 2? vi j c+ 2 tvi \ C \ 


1 



— h 


WpY\t (It 


1 


7T j- » (t 2 + - k*)(t 2 - &*)*(* 


i 


cop V 


1 


2 (z + k) 


TV 


Hz) — Hog) 

z — a i 


+ 


- z) 

Hz) 


2(2 + A) 

H — “0 


2 + «i 



Here <p(z) is given by the following expression, 


Hz) = (A 2 — 2 2 ) » arctan [(A - z)*(A + z)-»] ( 

while L+(z) is that multiplicative factor of L(z) which is regular in the upper 
half plane Im z > — A 2 and Z/ + (z) is its derivative, cti = (A 2 — a > 2 p 2 FJ)*, and 
it is the root which has a positive imaginary part. An equivalent form for 
Hz) which is useful for |z| -> °c, Im z appropriately limited, is 



For 


OC 


lin z > — k 2 we have that 




+ 0 (*- 2 ). 


We have finally upon integrating that 


£+(z) = (z + k)~t exp (— 


TV 



W) - *(*) - xP(~ ai ) 


t — Oil 


t + a i 



dt 


up to an arbitrary multiplicative constant. By a similar calculation in the 
upper half plane we find that L_(z) = L + (-z). The decomposition of M(z) 
is simdar to that of L(z) save for the fact that Y, is replaced by F 2 Hence 

K+(Z) ' S Simply the ratl ° of M Hz) to L + (z). We shall record here, for future 
use the asymptotic forms of K + (z) and K.(z) for z -> oo and Im z suitably 

limited. Using the asymptotic forms we gave for *(z), it is clear that 


K + (z) =c I 1 + ^ 


00 > Im z > — A' 2 , while 


- (Y 2 - FO — 

w z 


K-(z) = c, Tl - ^ (F 2 - F,) — 

L 7T z 


*’ Im 2 < The instants c x and c satisfy the condition that cc, = 


1. 
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6. The Fourier-transform Solution of Equation (13). Equation (13) may 

now be written in terms of K_(z) and K + (z) as 


(14) - 


3*i( 2 ) 1 + A i 1 1 

K+(z) t(z - k x ) iK^z) ~ K + {k x ) J 




where the left side is regular in the upper half plane Im 2 > -Im k, while the 
right side is regular in the lower half plane Im 2 < -Im 8 or Im k x . Since 
both sides are regular in a common strip, the left side is the analytic continua¬ 
tion of the right side and therefore each side of equation (14) is regular in the 
entire 2 -plane, i.e. y each side is an entire function E(z). The determination of 
E(z) in this case is a straightforward task since kJ(z) and K+{z) are bounded 
for \z\ > 00 , Im z suitably limited, and since 4>i(z) and 4>o(z) approach zero for 

\ z \ 00 > I m 2 suitably limited. Indeed from these facts we note that E(z) 

approaches zero for \z\ —> <», and an immediate application of Liouville’s 
theorem on bounded functions tells us that E(z) vanishes identically. This 
gives us then 



4> 0 (z) 


_ 1 Ai _ 

i(z — k x )K_{z)K+{k x ) 


From equation (15) we see that <f> 0 (z) = 0(zr l ) as \z\ —» 00 , Im 2 < Im k Xf 
and therefore </>o(x) = 0(1), x —> 0 + . Incidentally 4>o(z) is regular in the lower 
half plane Im 2 < Im k x , since the next singularity of $o(z) is k. If we 
examine the expression for $1 ( 2 ) for ( 2 ) —> 00 , Im 2 > —Im k, we find that 


#,(«) = - 1 . + A1 - J_ ± Ai 

u ’ iz izciK+{k x ) 
up to terms of 0(1/2). This tells us that 

«*)--(! +A0+±+±, (x —* O'). 

But 0i (x) = <t>(x,0) — (1 + A\)j so that for x —> 0“ we have that 



1 +iii 

ciK+(k x ) 


A similar calculation for <j> 0 (x) gives us the same result so that we see that 
0(x,0) is continuous at x = 0. 

7. The Far Field. The calculation of the far field, which incidentally will 
provide us with the diffraction pattern, does not offer any major difficulty. 
Care is required in the neighborhood of the angle 6 = — 0', and this brings us 
again to a problem which has occupied workers in diffraction theory since the 
days of Sommerfeld’s [6] initial contributions. It appears that if the complex 
integral is cast into an appropriate form, the problem of calculating the far 
field in the neighborhood of d = — O' is considerably simplified. 
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We start by considering representations of <fi(x,y) for x ^ 0. In the first 
place we have from the convolution theorem the following Fourier-integral 
representation for <t>(x,y): 


<t>(x,y) = <F(x,z/) + A\<t> r {x,y) 



T - A 1 ) 



K+ l {kz) exp [izx - iy(k- - z 2 )i] dz 
c K-(z)[icpY ! -f (k' 2 — z 2 )*](z — k x )’ 


where the path C is drawn in the strip of regularity — Im k < Im z < Im k z . 

Useful forms of this integral may be derived by deforming the path C into 

a hyperbolic path in the same fashion which we indicated in Sec. 3 . We write 

x — R cos 0, y = R sin 6, — 7 r < 6 < 0, R > 0, while z = k cos (0 -f- i T ), 

-oo < T < co on the hyperbolic path. The integral representation for 
<t>(x,y) then becomes 

<t>{x,y) = <F(x, 2 /) + D<f> (T) (x,y) + - --- - ^ 

Z7T 

f * _ exp (ikR cosh r)K^}(k x ) sin (0 + i T ) dr 

J- « K-[k cos (0 + irj][upYi — k sin (0 + zr)][cos (0 + zV) — cos 0'] ; 

where D = Ai for — tv < 0 < — 0' and D = A 2 for — 0' < 0 < 0. For 
0 = —d' the integral is defined in the sense of a Cauchy principal value. 

We now turn our attention to the integral 



_ exp (ikR cosh r)K+ l (k x ) sin (0 + ir) dr 

K~[k cos (0 + zV)][wp Y i — k sin (0 + zV)][cos (0 + ir) — cos $'] 


in the limit that \kR\ —» oo. 


This integral may be rewritten as 


I _ exp (ikR cosh r)K+ l (k x ) sin (0 - j T ) dr 

J - * cos (0 — it)][ cop Y I - k sin (0 — ir)][cos (0 — ir) — cos 0'] 

and in turn may be rewritten as the sum of the following two integrals: 


/ [/M +/( — r)] exp (ikR cosh r) 

( 16 ) / _ [upYijsc - sc/A) - ks 2 A(c - c'A ) - kcB 2 (A - cc')] dr 

J - ™ (A — cc' — ss')(A — cc' + ss')(o) 2 p 2 Yl ~ 2ksAwpYi -f- k 2 s 2 + k 2 B 2 ) 

+ [ * t/( T ) ~ /(~t)] exp (ikR cosh t)[g) P Y x ( A - cc ') + /c$(c 2 - A 2 )]B dr 
7 - * (.1 — cc' — ss')(A — cc' + ss')(u 2 p 2 Y 2 — 2ksAupYi + k 2 s 2 + k 2 B 2 ) * 

The following abbreviations have been employed in the last two integrals: 


/W = 


K?(k x ) 


K_[k cos (0 + zV)] ; 


5 — sinh r, c = cos 0, c' = cos 0' 


^4 = cosh t j 

s = sin 0, s' = sin 0'. 


We observe that each of the integrands in (16) is well behaved save possibly 
for the factor \/{A - cc' + ss'). For r = 0 and 6 = - 0 ’ this factor is singu- 
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lar. (The case 0 = 0' cannot occur here.) In order to isolate this difficulty, 
we rewrite (16) as 


1 / [/(t)+/(— t)]<7(t) - 2g(0)/(0) 

2 / _ „ .4 - cc' + ss' P ( tkR C0&h 7 ) dr 


+ 1 

9 



00 


2fi r (0)/(0) exp {ikR cosh t) 


— 00 


A — cc' + ss' 


dr 




00 


[/to -f(~r)]h(r) .. , , , 

^ ^ X P (^kR cosh r) dr 


where now 


and 


/ \ _ (sc — sc'.4)copFi — ks 2 A(c — c'A ) — A,*c£ 2 (^4 — cc') 
(A - cc' - ss')(«yF? - 2ksAu> P Y 1 + k 2 s 2 + k 2 B 2 ) 


h(j) = 


[c o P Y 1 (A - cc ') + ks(c 2 - A 2 )}B 


(.A - cc' - ss')(w 2 p 2 Y{ - 2ksAwpY 1 + k 2 s 2 + k 2 B 2 ) 


The integrands in the first and third integrals are continuous for r = 0, 
0 — — O', and an application of Kelvin’s method of stationary phase shows us 
that for 0 9 ^ —6' each of these integrals is 0(l/kR ) as \IcR\ —» oo. For 0 = —0' 
we get a contribution 0(l/y/kR) for \kR\ —> oc, but we shall not pursue this 
calculation further. 

This leaves us with the second integral, which may be rewritten as 


g(0),m 



OO 


[1 — cosh (r/2)] exp (ikR cosh r) dr 


— oo 


A — cc' + ss' 

+ g(0)f(0) 



00 


— 00 


cosh (r/2) exp (ikR cosh r) dr 

A — cc' + ss' 


The first of these integrals again provides us with a contribution 0(l/kR) as 
\kR\ —> co for 0 —O' and 0(l/\//W?) as |A*Z2| —> oc for 0 = —O'. The last 

integral may be expressed in terms of a Fresnel integral as 

— i \Ztt exp (i7T/4) exp [ikR cos (0 + 0')] 

|sin (0 + 0 ')/2\ 

This integral provides us with the basis of determining the behavior of 
<j>(x,y) as \kR\ —♦ oo. If 0 ^ —O', the dominant term is 0[exp (ikR)/y/kR] 
for \kR\ — > 00 . In the neighborhood of 0 = — O', we simply have a bounded 
contribution. Such a situation also arises in the transition from the illuminated 
to the shadow region in the excitation of an infinite half-plane wave in elec¬ 
tromagnetic or acoustic diffraction and has been discussed by Sommerfeld 
and others. In summary, we have that <f>(x,y) is asymptotic to 

4>(x,y) = <t>Kx,y) + Ai<l> r (x,y) 

exp (ztt/ 4)(/3 2 - j8 2 )y(0)/(0) V* exp [ikR cos (0 + 6 ') } [ " e <w ' dw 

27 r|sin (0 + 0 ')/ 2 | J\2[6in* (o+o')/2]kR\ 



oo 


e iwt dw. 


|2[sin 2 (6 + 0 ')/ 2 ) kR I 
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for |fc/2| —» °°, |0| > 0' while it is asymptotic to a similar expression for |0| < 0' 
with the constant A i replaced by A 2 . For 0 = — 0' other contributions arise, 
as we have already noted. 

Bibliography 

1. Albert E. Heins and Herman Feshbach, The coupling of two acoustical ducts , J. Math. 
Phys. vol. 26 (1947) pp. 143-155. 

2. Phillip M. Morse and Richard Bolt, Sound waves in rooms, Rev. Mod. Phys. vol. 16 
(1944). 

3. Raymond E. A. C. Paley and Norbert Wiener, The Fourier transform in the complex 
domain, Amer. Math. Soc. Colloquium Publications vol. 19 (1934) Chap. IV. 

4. Edmund T. Copson, On an integral equation arising in the theory of diffraction, Oxford 
Quart. Math. vol. 17 (1946) pp. 19-34. 

5. Edmund T. Copson, The asymptotic expansion of a function defined by a definite integral 
or contour integral, Admiralty Computing Service, 1946. 

6. Arnold Sommerfeld, Mathematische Theorie der Diffraction, Math. Ann. vol. 47 (1896) 
pp. 317-374. 

Carnegie Institute of Technology, 

Pittsburgh, Pa. 

Massachusetts Institute of Technology, 

Cambridge, Mass. 




ON THE DIFFUSION OF TIDES INTO PERMEABLE ROCK 


BY 

G. F. CARRIER AND W. H. MUNK 1 

1. Introduction. In the irrigation wells of the Hawaiian Islands, water- 
level fluctuations whose frequency components correspond to those of the 
ocean tides have been observed by Cox [1], In this paper we shall formulate 
the problem associated with this phenomenon, assuming the observed ground- 
water fluctuations to represent a diffusive transmission of the tidal disturbances 
through the porous volcanic structure of the island. The agreement with 
observation is fair. The results obtained here may find application in a 
method for estimating subsurface permeability. 


C 


x 


Fig. 1. Below BOC the medium is permeable. The pressure variations on OB depend only 
on the tidal water-level variations (mean water level: OA). 

2. The Formulation of the Problem. Consider the situation in Fig. 1. If 
the ocean level remained fixed, the pressure distribution in the fluid would be 

Po(y ) = -pogy. 

However, when the fluid is in motion, the pressure distribution depends on the 
velocity field. Within the porous medium, we adopt Darcy's law 




_ m _ 

Po 


grad (p - po), 


where p is the pressure, k the permeability of the medium (taken as isotropic 
and constant), p the viscosity of the fluid, m the mass flow per unit area, po 
the mean fluid density, and v as defined by m/p 0 is the “velocity.” 

Equation (2) replaces the law of conservation of momentum but we must also 
invoke the mass-conservation law. This takes the form 



div m 
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where A is that fraction of the volume which is filled with fluid (i.e., the 
porosity). We adopt a simple compressibility law for the rock and fluid, viz., 

(4) pA = p 0 A 0 [l + y (p — Po)], 

where y is essentially (p 0 c 2 ) -1 with c the speed of sound in the fluid. We may 
now combine the foregoing equations to obtain 



where q = p — p 0 and A is the Laplace operator. 

In equation (5) the expression k/pyA 0 has the dimensions of diffusivity 
(square centimeters per second), and the equation is of the diffusion type 
rather than a wave equation. The essential reason for this is that the pressure 
gradient is proportional to a velocity rather than to an acceleration. 

.The boundary conditions must now be given. On OB, the sea bottom, we 
may associate the “pressure” q with the tidal-wave amplitude. In fact, if 
we decide to consider a single periodic tidal-wave contribution, we write 


q(onOB) = q ie iut . 


On the free surface near Ox, we have from equation (2) 



where 7 j(x,t) is the ^-coordinate of the free surface. On this surface, however, 
p = 0, and using equation (1), 

00 q(x,v) = -po = pogv- 

We may now use equations (6) and (7) to obtain 



In these equations, q should be evaluated at y = rj but, in the usual linearized 
way, we apply equation (8) at y = 0. 

It is now convenient to define some dimensionless quantities, viz.: 


We then have 



pAcj y , x 

y ^k ~ L ~ y ’ L~ 

( _ plg 2 ky 
pA 06 ) 


(9) 

( 10 ) 
( 11 ) 


+ 9W = 
q(OB) = e", 

Qt(x',0,t) + gy (z',0,r) = 0. 
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If we now omit the prime on x and y and consider the periodic response of the 
system, we may write 

q(x,y,t) = <p(x,y)e iT 
and equations (9) to (11) become 

(12) — te<p = 0, 

(13) <p = 1 on OB, 

(14) <Pv(x, 0) + iip(x, 0) =0, (x > 0). 

This problem will be solved for two orientations of the line OB: (1) when the 
ocean-land interface is a vertical cliff (0 = 90°); (2) when OB is horizontal 
(0 = 0°). The latter case is the more realistic one since the bottom slope 
seldom exceeds 5 to 10°. 

3. The Horizontal Interface. In this problem, equation (13) takes the form 

(15) <p(x, 0) = 1, (x < 0). 

It is convenient to introduce the Fourier transform of <p and, in fact, we shall 
find (p by invoking the Wiener-Hopf technique. We define 

(15) £(£,*/) = f_^ <p(x,y)e~ i l x dx, 

and the transform of equation (12) becomes 

(17) Ipyy ~ (£“ + it) $ = 0. 

Hence 


( 18 ) jp = ^(Q^VF+E 

Now we can define hi(x) and h 2 (x) by 


h i(x) = lim 

a—* 0 


e° : 

0, 


and hi,(x) + h 2 (x) 



h,{x) = j 
However, 



(x < 0), 
(x > 0), 

(* < 0 ), 

(x > 0), 


Mi) = e(f,o) = ha) + a,($). 

It can be ascertained that 

(!9) £(1) = —U- 

a — 

and that h 2 (0) being the transform of a function which vanishes for x < 0, is 

regular in £ whenever Im (£) < 5, where 8 is a nonnegative real number. We 
can now define 

f( x ) = + i<p(x, 0 ) 


/(£) = (VC + u + *)ii(0 

= (Vc+Te + i)[hl(i) + 


and 

( 20 ) 
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Because of our boundary conditions, / (like A,(£)) is analytic in the upper half 
£-plane. We can anticipate, in fact, that each function of £ appearing in 
equation (20) will be regular in a horizontal strip lying below t(«/2)» and above 
the smaller of —i(e/ 2)5 and —ia. In this strip, then, equation (20) is a mean¬ 
ingful equation. We now anticipate that G(£) = \/£ 2 + ie + i can be written 

G(£) = <?+(£)<?_(& 

where G+(£) is analytic and has no zeros above £ = -i(e/ 2)* and where G_(£) 

has the same properties below £ = i(e/2)K We can write equation (20) in the 
form 



/(£) 

G+(£) 


- G_(-ta)ft,(£) = A,(£)G_(£) + [(?_(£) - G_(—fa)]^!(£). 


2 -plane 


- A 

- ; 

to » 






o 


Ar r > 


Fig. 2. Integration contour for equation 
(23). 


The left side is analytic in the upper halt 
plane, the right side in the lower, and, 
by the usual analytic-continuation 
arguments, we may associate each side 
with an entire function. Since, as 

S —► 00 , both Al({) —> 0 and 1 /(?+(£) —> 0, 
and since/(£) should be bounded at », 
this entire function must be identically 
zero. Thus, 



G-(-ia) 

. gm) 



The free-surface level 2 is then given by 






when we go to the limit a = 0. We cannot evaluate this, of course, until we 
find (?_(£). 

4. The Factorization of G(£). We shall use a well-known technique [2] to 
obtain (?_(£)/(?_ (0). We note that 

(23) In G(£) = In G + (£) + In G_(£) = ± f dz 

Zwi J r z - £ 


for all points in T, provided T = T i + r 2 is as indicated in Fig. 2. We may 
now define (? + to be that contribution of equation (23) associated with T i 
and (?_ to be that associated with r 2 . It is clear that the appropriate analytic 
continuations of these functions into their respective half planes of definition 
are automatic. We want to compute, then, 


i f+*+iVl72 

H(& = hi GM) = -«-./ In [G(z))(z - f )" 1 dz, 

Zwi j - oo+»vi72 


2 Corresponding to unit tidal amplitude. 
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and, for Im £ < 0, we may write 

(24) ff(£) = - g ZiJ_ m ln K 2 ' + «)* + *K 2 ~ 0" 1 dz ■ 

If we differentiate once with regard to £ and integrate by parts, we obtain a 
readily integrable formula for //'(£). In fact, the result is 


(25) H\z) = 


1 


+ c 


1 


2(z - iff) 2(z 2 + 0*) 


£ ln 


0 + 1 


ln 


2 + (z 2 + a 2 )* 


/i - 1 (z 2 + a 2 )* 2 - (z 2 + a 2 )* J 

where the branches must be chosen so that G'_ is regular in the lower half plane. 
In this formula a 2 = ie, (3 2 = 1 + ie. In the limit e —> 0, we have 


(26) 




+ lu 


i) 


-i 


1 ln f£ 

r e + 1 


The quantity GL(0)/GL(£) is defined by 


(27) 


G-(0) 

G-(£) 



( 


= - H’(z) dz, 


0 


where //' is given by equation (25) for the compressible situation 
equation (26) for the incompressible case. Equation (22) now yields 


and by 


(28) 


ht{x) =~ £ 

-7r / - * 


-i 




exp 


- / i H'(z) dz + i£x 
oj 


^ - 


1 


d£. 


The evaluation of this integral will be discussed in Sec. 6. 

6. The Vertical-cliff Problem. If we now consider the situation where the 
land-ocean interface is vertical, one must find the solution of 


(29) 

(30) 

(31) 


(in x > 0,y < 0). 

(y < 0). 

(x > 0 ). 


A<p = ie<p, 

<p(0,y) = 1, 

0) + i<p{x, 0) = 0, 

It is convenient to define \p(x,y) = + ip. The function i p must always 

obey equation (29), and the boundary conditions on \f/ become 

(32) 


H0,y) = i, 


and 

(33) 


C y < o), 


Hx,o) = o, 


(y > o). 


We can now extend the domain on which i p is defined to the entire right half 

plane and replace the present problem by that of finding +{x,y), which obeys 
(15) and the condition 

* = ieav < (y < 0), 

~* = J™ ■- ie ~ av > (y > 0). 


(34) 


i(0,y) = 
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Since will be odd in y, equation (33) will be automatically satisfied. The 
solution is now most readily obtained by using Fourier transforms again. We 


/_ „ dy. Multiplying equation (29) by e r^w and 


integrating yields 
so that 


hr — (y 2 + ie)$ = 0 
$ = A(y) exp (- vV + itx). 


The function A must be obtained from the boundary conditions, i.e., 

A = hO,y) = —2? 7(?? 2 + a 2 ) -1 . 

Finally, since $ = i(y -f- 1 )^, we have 


(35) 


1 


00 


* - s; - . - S /.. ^ ft £ * * +■ 


(v 2 + a 2 )( n + 1 ) 


In particular, on y = 0, we have 
(36) <t>(x,Q) = exp {—fix) — - 

7 T 



00 


ze 


xxz 


- . a/z 2 + a 2 (z 2 + a 2 + 1 ) 


dz 


= exp (-fix) + 


2 i 


7 T 



arcsinh - ) sinh fix 

a 


— I cosh [fi(x — t)]K 0 (<xt) dr 
o 


where K 0 is the modified Hankel function and fi 
limit e —> 0 , we note that 


= (1 + ie)K Finally, in the 


and 


arcsinh - 

a 


In ( 2 e-i) - 


ITT 

T 


a 


Thus, for 


K 0 {oi,r) ~ — In ^ — 7 — In r. 

0 , 

2 i f x 

<p{xfi) = e~ x d-/ (7 + In r) cosh (x — r) dr, 

n Jo 


where 7 is Euler’s constant. This function, (p{x, 0 ), behaves asymptotically 
like — 2i/irx. 

6 . The Results. We shall reduce the integral implied by equations (26) 
and (28) and those given by equation (36) to a form convenient for computa¬ 
tion. For the horizontal-interface problem with a = 0 (the incompressible case) 


(37) *(*,0) = 


2tt 



00 


[^exp |tt “ 1 ’ [In iz/{z 2 + 1 )] dz + i£ej — 1 J dz 


— 00 


= 2~^e~ x ~ ,T/8 T - 6ti/8 



Vi + Tz 

1 e~ zv — y~h~ x/v 

o (1 + y)K 1 - y) 

i __ 


COS M dz/ 



T • /» I e~ zv + y-h~ z/v . , 

- it Zg—tx/s / -r sin u dy, 

o (i + y)K l - y) 
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0 1 2 3 4 5 


X 

Fig. 3. Ground-water amplitude vs. distance inland. Curve I, incompressible horizontal 
interface. Curve II, incompressible vertical interface, x, e = 0.01 (vertical interface). 
O, average of observed amplitudes. K- 12, 12-hr. tide observation range at Kuau. K- 24, 
24-hr. tide observation range at Kuau. 



x 

Fig. 4. Ground-water phase lag vs. distance inland. Curve I, incompressible horizontal 
interface. Curve II, incompressible vertical interface, x, c = 0.01 (vertical interface). 
K-12, 12-hr. tide observation range at well. K- 24, 24-hr. tide observation range at well. 

where u = tt“ 1 [In v/(l — p 2 )] dv. The numerical integration of this latter 

form yields the results depicted in Figs. 3 and 4. 

1 he vertical-cliff result as given by equation (36) can also be written 



Asymptotically these functions [i.e., the *>(x,0) of equations (37) and (38)] 
behave like -i/tx and respectively, for a = 0. The results of the 
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numerical integration of this equation are also plotted in Figs. 3 and 4, both 
for 6 = 0 and e = 0.01. The notation on these curves is given by 

(39) <p(x,0) = B(x)e~ itl(x) . 

The experimental observations which are available are scanty, and the spread 
of the data is great. The observations at the Kuau well [1] are also indicated 
in Figs. 3 and 4. Note that the phase-lag observations are not easily inter¬ 
preted but that the amplitude decay is reasonably consistent. Since, for a 
given well the “distance” x is proportional to the frequency, the 12- and 
24-hr. tides should correspond to x values having the ratio 2. Choosing the 
location of these data points on this basis, we find that this location is consistent 
with a permeability value of 5 X 10' 6 cm. 2 and a value of € of order 0.01. This 
permeability is consistent with other estimated values for this volcanic mate¬ 
rial. The observed phase lag,'-however, is much too low to be consistent 
with our theory. Furthermore, as is indicated by the vertical-cliff computa¬ 
tions, the compressibility correction would lead to a greater discrepancy—not 
a smaller one. Thus, the phase-lag results are somewhat disappointing 
although the amplitude decay is consistent with our theory. This is not too 
surprising since the nonhomogeneity of the geologic setup can easily upset the 
validity of the results. Nevertheless, it is difficult to explain the very low 
phase-lag observations and especially the fact that this phase lag seems rela¬ 
tively independent of frequency. 

Finally, it seems of some interest to record the asymptotic behavior of the 
function defined by equation (28) for the compressible case. One can show 
that 

A*(») ** - £ Ki(ax) - J" K,(t) dr - [K\{ a x) + K 0 (<xx)] + • • • , 

where K y is the Hankel function of imaginary argument. However, this series 
is an asymptotic representation in 1 /a rather than in x, and the function is well 
represented only when x » |ln a\. That is to say, it is useful only in a range 
of x where the compressibility is important. Since \a\ = 0(0.1) and x = 0(1) 
in the cases of interest to us, these results are not particularly pertinent to the 
present problem. 
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SOME REMARKS ON RADIATION CONDITIONS 


BY 


J. J. STOKER 


In linear wave-propagation problems for the so-called steady state (to be 
understood here as a motion that is simple harmonic in the time) in unbounded 
domains it is in general not possible to characterize uniquely the solutions 
having the desired physical characteristics by imposing only boundedness con¬ 
ditions at infinity. It is, in fact, necessary to impose sharper conditions. In 
the simplest case in which the medium is such as to include a full neighborhood 
of the point at infinity that is in addition made up of homogeneous matter, the 
correct radiation condition is not difficult to guess. It is simply that the wave 
at infinity behaves like an outgoing spherical wave from an oscillatory point 
source, and such a condition is what is commonly called the radiation, or 
Sommerfeld, condition. Among other things this condition precludes the 
possibility that there might be an incoming wave generated at infinity— 
which, it not ruled out, would make a unique solution of the steady-state 
problem impossible. 

If the refracting or reflecting obstacles to the propagation of waves happen 

to extend to infinity—for example, if a rigid reflecting wall should happen to 

go to infinity—it is by no means clear a priori what conditions should be 

imposed at infinity in order to ensure the uniqueness of a steady-state solution 

having appropriate properties otherwise. 1 The point of view to be presented 

here is that the difficulty arises because the steady-state problem is an unnatural 

problem in mechanics and that, in principle at least , one should rather formulate 

and solve an appropriate initial-value problem and then find the solution of the 

steady-state problem by making a passage to the limit in allowing the time to tend 
to infinity . 2 

The steady-state problem is unnatural because one makes a hypothesis about 
the motion that holds for all time, while Newtonian mechanics is basically con¬ 
cerned with the -prediction —in a unique way, furthermore—of the motion of a 
mechanical system from given initial conditions. Of course, in mechanics of 
continua that are unbounded it is necessary to impose conditions at infinity 
(in the space variables, that is) not derivable directly from Newton’s laws, but 

for the initial-value problem it should suffice in general to impose only bounded¬ 
ness conditions at infinity. 


1 For a treatment of the radiation condition in such cases see F. Rellich, Uber das asymp- 

tohsche Verhalten der Losungen von Au + ti =0 in unendlichen Gebieten, Jber. Deutschen. 

Math. Verein vol. 53 (1943). Also see F. John, On the motion of floating bodies II, Comm 
Pure Appl. Math. vol. 3 (1950) p. 54. 

2 Ihe formulation of the usual radiation condition is doubtlessly motivated by an instinc¬ 
tive consideration of the same sort of hypothesis combined with the feeling that a homo¬ 
geneous medium at infinity will have no power to reflect anything back to the finite region. 
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If one wished to be daring one might, on the basis of these remarks, formulate 
the following general method of obtaining the appropriate radiation condition 
for the steady-state problem: Consider any convenient problem in which the 
part of the domain outside a large sphere is maintained intact and initially at 
rest. (In other words, one might in particular feel free to modify in any con¬ 
venient way any bounded part of the medium.) Next solve the initial-value 
problem for an oscillatory point source placed at any convenient point. A 
passage to the limit should then be made in allowing the time t to approach 
infinity, and after that the space variables should be allowed to approach 
infinity. The behavior at the far-distant portions of the domain should then 
furnish the appropriate radiation conditions independent of the constitution of 
the finite part of the domain. It might be worth pointing out specifically that 
this is a case in which the order of the two limit processes cannot be inter¬ 
changed: obviously, if the time t were first held fixed while a space variable 
tends to infinity, the result would be that the motion would vanish at infinity* 

because of the finite propagation speed for disturbances, and no radiation con¬ 
ditions could be obtained. 

The writer would not have set down these remarks—which are of a character 
so obvious that they must also have occurred to many others—if it were not for 
two considerations. Every reader will doubtlessly have said to himself, “That 
is all very well in principle, but will it not be prohibitively difficult to carry out 
the solution of the initial-value problem and to make the subsequent passages 
to the limit?” In general, such misgivings are probably all too well founded. 
However, the writer has come upon an interesting special case in which (1) the 
indicated program can be carried out in an elementary way in all detail and (2) 
it is slightly easier to obtain an integral representation for the solution of the 
initial-value problem than it is to solve the steady-state problem with the 
Sommerfeld condition imposed. 

The problem in question concerns the propagation of surface-gravity waves 
in an infinite ocean. We consider only irrotational incompressible flow and 
assume further that the motion is a small oscillation in the neighborhood of the 
rest position of equilibrium; a linearized theory can then be obtained. We 
restrict ourselves to two-dimensional motion in an xy- plane, with the y-axis 
taken vertically upward and the .r-axis in the originally undisturbed horizontal 
free surface. Under these assumptions there exists a velocity potential 
<t>{x,y,t) which is a harmonic function in the lower half plane: 

(1) 0xx + 0 yy = 0, (y < 0 ,t > 0). 

(Here and in what follows subscripts denote partial derivatives.) The free- 
surface boundary conditions (cf. Lamb, Hydrodynamics, p. 364) are 

(2) — 0„ + Vt = 0, 

(3) 4> t + gv = - - Pi 

p 

In these equations rj = rj(x,t) represents the vertical displacement of the free 
surface measured from the rc-axis, and p = p(x,t) represents the pressure applied 
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on the free surface. We suppose that <£ and its first and second derivatives 
tend to zero at «> for any given time t —in fact that they tend to zero in such a 
way that Fourier transforms exist—but we do not, in accordance with our dis¬ 
cussion above, make any assumptions about the behavior of our functions as 

At the time t = 0 we prescribe the following initial conditions, 

(4) <£(£,0,0) = <£*(£,0,0) = p(x,0) = 0, 


which imply that the free surface is initially at rest in its horizontal equi¬ 
librium position. 

In what follows we consider only the special case in which the surface pres¬ 
sure p(x,t) is given by 


(5) p(z,t) = 5(x)e iut , (t > 0). 

Here 5(x) is the Dirac 5-function; in other words, an oscillatory pressure point 
is applied at the origin on the free surface and maintained there for all time. 
By inserting this expression for p in (3) and eliminating the quantity rj by 
making ue of (2) the free-surface condition is obtained in the forms 


(6) g<t>u + <t>u = — — 8(x)e iut , (t > 0). 

Our problem now consists in finding a solution <t>(x,y f t) of (1) which behaves 

properly at <*> and which satisfies the free-surface condition (6) and the initial 
conditions (4). 3 

We proceed to solve the initial-value problem by making use of the Fourier 
transform applied to the variable £, or rather, since the problem is evidently 
symmetrical with respect to the y- axis, we may use the cosine transform. 4 * 
The result of transforming (1) is 

(7) — s 2 <£ + <£„„ = 0, 

in which 4>(s,y,t) is the cosine transform of <t>(x,y,t) and use has been made of 
the fact that <f> x = 0 for x = 0, y < 0, and <£ x —> 0 as x —■> <x>. The bounded 
solutions of (7) for y < 0, s > 0 are all of the form 

( 8 ) 4>(s,y,t) = A(s,t)e 8y . 

The cosine transform is now applied to the boundary condition (6), with the 
result 


3 The steady-state problem is formulated as follows: One sets <f> = \p(x,y)e iut . The time 
factor e iut can be canceled from (1) and (6), and the harmonic function then satisfies the 

free-surface condition g\p v — a, 2 ^ = — —- 6(x). In addition one would postulate that the 

motion should have the character of an outgoing progressing wave at x = — oo and x = + ». 

This problem was solved long ago by H. Lamb, On deep water waves , Proc. London Math. 

Soc. Ser. 2 vol. 2 (1904) pp. 371-400 (see pp. 388-389). It is a curious fact that this problem, 

which involves the determination of a harmonic function and not a solution of a wave 

equation, nevertheless has not a unique solution of the desired type if only boundedness 
conditions are imposed at 

4 The same general approach to problems of this type is used in the book by I. N. Sneddon, 

Fourier transforms, McGraw-Hill Book Company, Inc., New York, 1951, p. 278. 
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9<i> v + <t>u =- \= — e iat , 

V2 r P 

and on substitution of 4>(s,0,t) from (8) we find 

Att + gsA =-— e iut . 

\/2ir P 

The initial conditions (4) now furnish for A(s,t) the conditions 

(11 ) 4(s,0) = 4,(*,0) = 0. 

The solution of (10) subject to the initial conditions (11) is 


(for y = 0), 


( 12 ) 


^l(s,f) = - 


1 id) 


y/%t P 


t p\ 0){t—T) 


o y/gs 


y/gs 


Finally, we insert the last expression for A(s,t) in (8) and apply the inverse 

transform to obtain the following integral representation for our solution 
<t>(x,y,t): 


(13) 


4>(x,y,t) = - 


l CO 
P7T 


» 


e sy cos sx 


t pi(i)(t—T) 


0 


/-=- sin \/gs r dr ds 
o y/gs 


This solution is also seen to be the unique solution of our problem. Our 
object is to study the behavior of this solution as t —> oo. 

Since y is negative (we do not discuss here the limit as y 0, that is, the 
behavior on the free surface), the integral with respect to s converges well and 
there is no singularity on the positive real axis of the complex 5 -plane. How- 
evei, the passage to the limit t > oc is more readily carried out by writing the 
solution in a different form in which a singularity—a pole, in fact—then 
appears on the real axis of the 5-plane. (It seems, indeed, likely that such an 
occurrence would be the rule in any considerations of the present kind since the 
limit function as t > co would not usually be a function having a Fourier 
transform, and one could expect that the limit function would somehow 
appear as a contribution in the form of a residue at a pole.) It is convenient— 
and of course legitimate by Cauchy’s theorem—to replace the path of integra¬ 
tion along the 5-axis in the 5-plane by a path L indicated in the accompanying 
figure. Ihe path L lies on the positive real axis except for a semicircle in the 
upper half plane centered at the point 5 = w 2 /#. 
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We now replace sin V^s x in (13) by exponentials and carry out the integra¬ 
tion on x to obtain 6 


(14) <t>(x,y,t) = - 


icoe 


iut 


Tp 



e* v cos sx 


1 «)< 


l (> — i(y/os+u>)t 


2 Vas Vgs - co 2\/ss Vgs + w 


+ 


——il rfs - 

firs - cu 2 J 


We wish now to consider the three items in the bracket separately, and as we 
see, two of them do indeed have a singularity at s = co 2 /g which is bypassed 
through our choice of the path L. The first two items are rather obviously the 
result of the initial conditions and hence could be expected to provide tran¬ 
sients which die out as / —> qo . This is in fact the case, as can be seen easily 
in the following way: That branch of \/.s is taken which is positive on the posi¬ 
tive real axis, and we operate always in the right half plane. If, in addition, .s 
is in the upper half plane, it follows that i(y/gs ± co) has its real part negative 
(w being real). Consider now the contribution furnished by the first item in the 
square brackets. Since the exponential has a negative real part on the semi¬ 
circular portion of the path L, it is clear that as t —> + oo, this part of the path 
makes a contribution that tends to zero as t —» oc . The remaining portions of 
L, which lie on the real axis, are then readily seen to make contributions which 


die out like 1 ft: this can be seen easily by integration by parts, for example, or 
by application of known results about Fourier transforms. The second item 
in the square brackets has no singularity on the real axis, so that the path L 
can be taken entirely on the real axis; thus, in accordance with the remarks 
just made concerning the similar situation for the first item, it is clear that this 
contribution also dies out like \/t. Thus for large t we obtain the following 
asymptotic representation for </>: 


/1 r \ ,/ .\ 10) . t f e* y cos sx , 

( lo ) 4>(x f y,t)d' -e‘ wt / - -ds. 

Trp Jl gs — ar 

Aetually, the right-hand side is the solution of the steady-state problem—as 
obtained, for example, in the paper of Lamb cited previously—when the con¬ 
dition at oc is the radiation condition stating that </> behaves like an outgoing 
progressing wave. The steady-state solution furnished by (15) is a little more 
awkward to obtain directly through use of the radiation condition—at least 
by methods known to the author—than obtaining solution (13) of the initial- 
value problem. In particular, the asymptotic behavior of an integral repre¬ 
sentation must be investigated in this case also before the radiation condition 
can be used. Thus we have verified in this special case that the radiation 
condition can be replaced by boundedness conditions (in the space variables, 

that is) if one treats an appropriate initial-value problem instead of the steady- 
state problem. 

6 The argument which follows was suggested to the writer by his colleague A. S. Peters. 
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Even though not strictly necessary—since (15) is known to furnish the desired 
steady-state solution—it is perhaps of interest to show directly that the right- 
hand side of (15) has the behavior one expects for an outgoing progressing 
wave when x -> + °o. The procedure is the same as that used above in dis¬ 
cussing (14): The factor cos sx is replaced by exponentials to obtain 










By the same argument as above one sees that the first integral makes a con¬ 
tribution that tends to zero as x —» -j- <». The second integral is treated by 
deforming the path L over the pole s = tf/g into a path M which consists 
of the positive real axis except for a semicircle in the lower half plane. The 
contribution of the second integral then consists of the residue at the pole plus 
the integral over the path M. But the contribution of the latter integral is, 
once more, seen to tend to zero as a; —> -f~ 00 because of the factor Thus 

behaves for large x as follows: 


(17) H(o»*/Q)x—ut) 

9P 

This represents a progressing wave in the positive ^-direction which, in addi¬ 
tion, has the wavelength 2wg/ co 2 appropriate to a progressing sine wave with the 
frequency in water of infinite depth. 

New York University, 

New York, N.Y. 



ON THE THEORY OF SCATTERING OF PLANE WAVES 

BY SOFT OBSTACLES 1 


BY 

E. W. MONTROLL AND J. M. GREENBERG 

1. Introduction. One of the most important physical processes under 
investigation today is the scattering of waves by various types of obstacles. 
In the realm of basic physics and physical chemistry, experimental scattering 
techniques are used to study the detailed character of atomic nuclei (and 
nuclear forces), the arrangement of atoms in simple molecules, and the sizes 
and shapes of large molecules such as the virus proteins. Observations of the 
scattering of light from distant stars are used to deduce the nature of inter¬ 
stellar dust. 

The general procedure is to choose a model for the scatterer and to compare 
observed scattering patterns with those calculated on the basis of the model. 
If an observed pattern is in agreement with that deduced from the model, one 
has considerable evidence that it is a reasonable one. This approach is of 
course quite old, having been used by Lord Rayleigh in his theory of the color 
of the sky and by Rutherford in his theory of atomic structure. 

Modern technology is also concerned with scattering processes. For 
example, the scattering of sound and radio waves is important in the fields of 
long-range communication and the application of radar and sonar. 

Although scattering processes have been known and investigated for many 
years, the theory of the processes has usually lagged somewhat behind the 
demands made upon it. Even though a tremendous amount of work has been 
done recently to rectify this situation, there is still much to be done by applied 
mathematicians and physicists in this field. 

The aim of this paper is to report on some progress that has been made by 
van de Hulst, Hart, Glauber, and the authors on the problem of scattering of 
plane waves by “soft obstacles.” An obstacle is considered to be soft if the 
wavelength of the wave inside the scatterer does not differ much from that of 
the incident wave in the absence of the scatterer. The well-known approxima¬ 
tion of Rayleigh, Cans, and Born in which the internal field of an obstacle is 
taken to be that which would have existed in the absence of the scatter is 
ordinarily used in the analysis of scattering by soft obstacles. Unfortunately, 
when the range of the scatterer is more than a few wavelengths of the incident 
radiation, this approximation becomes very poor. 

2. Fundamental Equations. We shall introduce the fundamental equations 
needed for the development of scattering theory in this section. Although we 
restrict ourselves to scattering of scalar waves with “quantum-mechanical” 
boundary conditions—wave functions and their first derivatives continuous— 

1 This work was supported by the Office of Naval Research. 

103 



104 


E. W. MONTROLL AND J. M. GREENBERG 

one can (not always without some effort) generalize the methods and results to 

vector waves and other boundary conditions. These boundary conditions 

correspond to the acoustical case in which the density of the scatterer is the 

same as that of the medium. They are the scalar analogue of the scattering of 
light by a dielectric. 

The fundamental wave equation is 

VY + k\r)i = 0. 

The quantity |*|* has several physical interpretations depending on the 
problem at hand. It is the particle density in the quantum theory, propor- 

tional to sound intensity in acoustics, and to light intensity in the corresponding 
vector equation of light scattering. 

We let 


( 2 ) 


lAo = exp Hc 0 (t • s 0 ) 


represent the wave function of an incident plane wave propagated in 
direction of the unit vector s 0 . Its wavelength X 0 is related to k„ through 


the 


(3) 


k 0 = 


2tt 


In the quantum-mechanical case we have 


(за) 

( зб ) 



2 mE 


2 m[E - V(r)) 
h 2 


where E is the energy associated with an incident particle, m its mass, V(r) the 

potential of the scattering field, and h = Planck's constant/^. We shall 
frequently refer to 

(3c) m(r) = ~1 

A'o 

as the relative index of refraction of scatterer to medium at the point r and to 

rrt = ki/ko as the relative index of a uniform homogeneous scatterer to its 
medium. 

If Mr) is the wave function of the scattered field, the solution of (1) is the 
form 


(4) i£(r) = M r ) + exp z7c 0 (r • s 0 ) 

where \p s (r) —> 0 as 1/r. Hence \p 8 (r) satisfies the differential equation 

(5) VY. + Kh = [k 2 - k 2 (r)][f/ 

and therefore the equivalent integral equation [1] 

(6) UR) = ^ / CXP |R l R r |~ f| f* 2 (r) - imr) dr; 

the integration extends over all space. 
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A physical interpretation of ( 6 ) based on considering every infinitesimal 
element of the scatterer as a Rayleigh scatterer is suggestive in seeking an 
approximate solution of (6). The quantity (exp ik 0 \R — r|)/4rr|R — r| is 
proportional to the wave function of the radiation emitted from an infinitesimal 
sphere at the point r. The term fc 2 (r) — h'l expresses the difference between 
the scatterer and the medium at the point r, and \f/(r) is the internal field at the 
point r which excites the infinitesimal scatterer of volume dr. Hence ( 6 ) 
corresponds to a superposition of the scattered radiation of a large number of 
infinitesimal Rayleigh scatterers, each excited by the proper internal field at its 
point of location. Hence if on some physical basis we can make a good 
approximation to the internal wave function \p(r), we can substitute it into 
(fi), integrate, and obtain an approximate scattered wave function. 

The theory of scattering most often used, the Rayleigh [2]-Gans [3]-Born 
[4] approximation, is obtained by approximating \p(r) by the wave function 
associated with the incident wave in the absence of a scatterer. If in the 


quantum-mechanical case \E\ » |F(r)| so that |/c 2 (r) — k 2 \ is small, this approxi¬ 
mation is good unless the scattering center has a range larger than several 
wavelengths of the incident wave. When the range is large, the incident wave 
becomes more and more distorted with increasing distance of passage through 
the scattering field. An improved estimation of the internal field will be 
discussed in the remainder of this article. 

When observations are made at great distances from the scatterer, (fi) can 
be simplified to 



UR) 


gikoK I 

4 ^ / 'l'(r)[k 2 (r) - A* 2 ] exp (-z7c 0 r • Si) dr 


where Si is a unit vector in the direction of the scattered wave. Although 
improved scattered wave functions can be obtained in principle by substituting 
an approximate internal wave function into ((>) and iterating several times, it 
becomes an almost hopeless procedure in practice because \f/ 8 (R) can seldom be 


computed for all R > 0 by ( 6 ). The asymptotic \p s (R) as R —> oo is usually 
quite easy to derive from (7), but unfortunately even a second iteration for 
UR) requires yp(R) for all R. 

Since the ingenious variational principal of Schwinger and Levine is essen¬ 


tially a second iteration process, it is difficult to apply in most cases. 


3. Scattering by Soft Spheres. We consider a sphere to be soft when the 
internal propagation constant hi (h\ = k (r) for r < a, the sphere radius) is 
almost equal to the external constant k 0 . In quantum-mechanical language a 
spherical obstacle corresponds to a square-well (or barrier) scattering potential. 
The soft sphere represents a well or barrier whose depth or height is small 
compared with the energy of the incident wave. The R.G.B. approximation 
applies very well to the scattering by a soft sphere whose radius is not too much 
laigei than the wavelength of the incident wave. When several wavelengths 
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of the incident wave fit into the sphere the incident wave function is a bad 
approximation to the sphere’s internal wave function. 

The qualitative nature of the internal wave function of the sphere can be 
guessed by considering the analogous one-dimensional problem of the propaga¬ 
tion of a plane wave through an infinite flat plate. Let the plane wave be 
propagated from left to right in Fig. 1. The wave will be partially reflected 
by the surface of the plate. The transmitted wave will also be partially 
reflected by and partially transmitted through the second surface of the plate. 

Hence, the forms of the wave function for the three 
regions of Fig. 1 are: 

I. i p(x) = e ixk ° + Ae~ ixk °, 

II. \p(x) = Be ixki + Ce ~ ixkl , 

III. \P(x) = De ixk \ 

where the constants A, B, C, and D are to be determined 
from the boundary conditions \p(x) and yp'{x) continuous, 
at the boundaries of the plate. 

One might expect the internal wave function of a soft sphere to be of the 
form of a linear combination of a transmitted and reflected wave with the 
propagation constant k\ associated with the scatterer. This is indeed the 
case, as we shall show below. 

The exact expressions for the internal and scattered wave functions of a 
spherical obstacle are easily derived from the fundamental wave equation 

2tt\ 

/ 

= 

~ Xo / 

(X 0 = wavelength of incident wave and Ai that of the corresponding wave in a 
medium composed of the material of the sphere) with the boundary conditions 




Fig. 1 


d(r\p) > (continuous at r = a). 


Let the incident wave be propagated in the 2 -direction. Then its wave 
function is 



= gizfCQ = — Q\k$r cos 0 


It is easily shown in the usual manner [5-7] that the transmitted and scattered 
waves are given by 

oo 

(10a) ^ A,,„i n (2n + l)j n (rki)P n (cos 6), 

n “0 

00 

(106) i. = ^ A,,J n (2n + l)h n {rko)P n (cos 9), 

n = 0 



where 

(lla) 

(116) 
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kih n (^ako)j n-+-i(^^i) koj i) hn+i(ciko) 

4 _ kij n (ak 0 )jn+i(aki) — k 0 jn(aki)j n+l (ako) 

kih n (ako)jn+i(.aki) — kojn(,aki)h n+ i(ak 0 ) 


Then 




+ i P,, 

(if r > a), 


i = 'Pt, 

(if r < a). 


Here z n (u ) is the spherical Bessel function 

2 n(«) = Q0 ^n+j(u), 

where z and Z represents either j and ./ (Bessel) or h and H (Hankel) function. 

The rapidity of convergence of (10a) and (101;) depends on the value of the 
ratio x = 2ira/\. Debye [8] has shown that the number of significant terms 
in these series is approximately the first integer larger than x. This result can 
be observed empirically in scattering tables [9, 10]. The following argument 
gives a physical basis to Debye’s theorem: If the series (10a) and (106) are 
interpreted as the internal and scattered wave functions of a particle scattered 
by a force center of range a, each term in the series represents an angular 
momentum state of the particle. The Ith state has an angular momentum 
h[l{l + l)]h Classically those particles which are outside the range a of the 
scattering force are unaffected by it. The angular momentum with respect 
to the force center of a particle of linear momentum p at the boundary of the 
range of the sphere of interaction is pa. Hence, the largest l which can con¬ 
tribute significantly to the series (10) is defined by pa = h[l(l + l)]b But a 
particle of momentum p corresponds to a wave of wavelength X = 2irh/p. 

Hence, / is related to X and a by 2 tt a/X = [/(/ + 1)]*, which for large l is Debye’s 
result. 

The Debye theorem suggests an approach to the approximation of the sums 

\ps and it. Let us suppose we can accurately approximate (11a) and (116) by 

simpler expressions when n < x so that the series for i, and i t can be summed 

with the approximate coefficients. If the approximate coefficients for n > x 

are small, even though they might be in considerable error, Debye’s theorem 

would indicate that the new summable series should be good representations of 
i, and \p t . 

We shall now attempt to simplify the coefficients (11a) and (116) by applying 
the asymptotic formulas [11] 

h n (x) ^ ( —i) n+1 x~ 1 e ix , 
j n (x) ^ x _1 cos -- —l l 


(12a) 

( 126 ) 
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which are valid when x » n. These formulas yield the following approxima 
tion to the denominators of (11a) and (116): 



+ ho) _ ia(k ._. 


2a 2 kik 0 


— e 


l ~ A ' #) [l + ( — ) n ke 2iaki ]; 



k\ + ko/ 


An especially interesting feature of this formula is that it reduces to the exact 
formula for the denominator as k x -» k 0 for all integral n’s. Hence, a result 
based on (13) should be valid for soft spheres, and indeed it should be exact in 
the limit as ki —» k 0 . Unfortunately, (13) is a very poor approximation to the 
denominators of (11a) and (116) in the case of |A| large. Then it yields coeffi¬ 
cients to (10) which do not satisfy the Debye condition of being small when 
n > 27ra/A. 

When the approximate coefficients, obtained by employing (13) in (10a) and 

(106), are used, both series can be summed by using well-known Bessel-function 

formulas. We first sum (10a) to obtain the internal wave function. For 
\r\ < a we have 


(14a) 




gia(ki—k 0 ) 

1 (A, + A*o)(l - **««■*>) 


5 > 


- ( — )"ke 2iak ']i n (2n + l)j n (rki)P n (cos 0). 


If we neglect the term proportional to A*, the sum is exactly the spherical 
harmonic expansion of a plane wave propagated in the 2 -direction, 


e ik ' lZ = exp {ik\r cos 0). 

I he term proportional to A: in the sum is the corresponding expansion for a plane 
wave propagated in the - 2 -direction. Hence 


(146) Ur) 


2A'iC' tt(A ' ,_Ao) 


(A'i + A*o)(l — k 2 e Aiakl ) 


exp [t(r • s 0 )A'i] 


ke liakx exp [ — i(r • s 0 )A’i] I 


so that the internal wave function corresponds to the superposition of two 
plane waves, one propagated in the direction of the incident wave and the 
other in the opposite direction. The vector s 0 is a unit vector in the direction 
of the incident wave. 

The scattered wave function can be computed with the appropriate denomi¬ 
nator in A KfTl by using the addition theorem for Bessel functions. However, as 
a preview of Sec. 5 we find \f/ 8 (r) by substituting (146) into the integral equation 
(7) with 



A'o 2 , 


(r < a), 

(r > a), 


k 2 (r) - k 2 = 


y 
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and, proceeding in the same manner as is done in the Born approximation, we 
obtain for large R 


(15) 

where 

(15a) 

and 

(156) 


) 2irR(l - AV“*') (/l ” ke 


h = J v exp [ir • ( k x s n — A 0 Si)] dr 
1 1 exp [ — ?'r • (AiSo + A 0 s,)] dr. 


The integration extends over the volume of the obstacle. 

To integrate 1 1 we define y as the angle between the fixed vector (/ciS 0 - 
and the variable vector r which ranges through the entire sphere. Let 


koSi) 


(16a) 


w i = 


/‘ i s o — A' 0 Si| 2 = k'l + k'l — 2/ci/vo cos 0, 


where 6 is the angle between s 0 and s Y . Then 


(166) 


/1 “ -tt exp (zrcoi cos y)r 2 sin 7 ^7 dr 

= 4 7 r 3 w sin a 


= 47ra 3 - 



du 


(a^i) Vj(a«i), 


where J%(x) is the Bessel function of order Similarly we have 


(16c) 

where 

(16d) 


/ 2 = 47ra 3 ( ^ ) (avi)~Ui(avi), 


»? = ** + A 3 - 2A,Ao cos 0. 


B.V substituting (166) and (16c) into (15) we can easily derive an approximate 

expression for the differential-scattering cross section a,(d) which has the 

property that a,(6) dii is the fraction of incident energy associated with the 

incident plane wave which is scattered through a solid angle dil at 6 We 
find [7] 


(17) v .{6) = R*U.\* 


= I A 


with 


«/j 2 (flq?i) 

_ (aon) 3 


A 


8 


- 2k cos (2aAi) •! A au d J i( av d . k 2 

(a 2 Vi wi)^ (aiu ) 3 

_ 47tA'J(A’i — /co) 2 a 6 


1 + k 4 — 2k 2 cos 4aA i 


Hus result is exactly that obtained by summing (106) with the approximate 
denominator (13) used in (116). 
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A similar result can be derived for the scattering cross section of an infinite 

circular cylinder. In the case of normal incidence of the plane wave to the 
cylinder one has [ 12 ] 


<r,(0) 


ns*/ 


2ira 4 ki(m — l) 2 fcg 
1 + k 4 + 2 k 2 cos 4 aki 


+ 2k sin 2 ah 


(acoi) 


a^iOi 



where on and V\ are defined as before with 0 as the polar angle about the 
cylinder in cylindrical coordinates. 

A somewhat different analysis of the problem of scattering by soft spheres 
has been made by van de Hulst [13] (see especially pages 45 to 46). First he 
applies Huygens' principle to the problem. In the limit as m = ki/k 0 -> 1 a 
negligible part of the incident radiation is reflected, and the refraction of the 
transmitted part causes only slight deviation in the direction of an incident 
lay. Hence, in quantum-mechanical language, a particle that passes through 
the sphere changes neither its direction nor intensity. It, however, changes 
its phase by an amount proportional to the length of the chord along which it 
passes through the soft sphere. As will be discussed at greater length in Sec. 

6 , the change in phase along a chord whose closest 
distance to the center of the sphere is u is 

x(u) = 2k 0 {m — l)(a 2 — u 2 )* (u < a) 

= 0 , (u > a). 

The wavefront for particles which traverse the 
sphere lag behind those which do not. The entire 
wavefront then develops a hemispheroidal depres¬ 
sion as it passes through the sphere. By a direct 
Fig. 2 application of Huygens' and Babinet's principles to 

the further development of the wavefront, van de 
Hulst] shows that if the scattered wave function is written as ^,( 0 ) = R~ ] f(6) 
exp (ik 0 R), then 



(18a) /(<?) = —a 2 ko (1 


e 2 ik 0 (m Doco» y)J 0 (z sin 7 ) sin 7 cos 7 dy, 


where u = a sin 7 and 7 is the angle between a chord and its associates normal 
to the sphere. The parameter z is defined by z = bk 0 sin 0 . 

It is interesting that this expression, which was derived by an application of 
Huygens’ principle, is equivalent to the limiting result as m —» 1 that follows 
from the exact formula (105) by applying the Bessel-function approximations 
(12a) and (125) to both the numerator and denominator of (105). On this basis 
van de Hulst showed that 


(186) f(0) = -aHk 0 | o ' /2 (l 


e -2(oMm-i)co5 7) f / 0 ( 2 sin 7 ) sin 7 cos 7 dy t 
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where now 



The difference of the phase factors of the two forms of /(0) comes from a differ¬ 
ence in definition of phase in the two analyses. In the limit as m —* 1 (18a) 
yields an expression for the differential-scattering cross section which is equiva¬ 
lent to the corresponding limit of (17). 

Equation (18a) can be written in a slightly different form if we choose the 
distance u of the chord from the center of the sphere as the variable of integra¬ 
tion and express the integrand in terms of x(u)> the phase changes at n. Then 

(186) f(6) = b f (1 - rx(“»)J o (koU0)u du. 

1 Jo 

4. The Approximating Sphere. Let us now consider a class of spherically 

symmetrical soft scatterers with the properties: 



ii. |< 7 (x)| is a monotone nonincreasing function of x which approaches 0 


sufficiently rapidly as x —> « so that xg{x) dx and x 2 g 2 (x) dx exist. 

The sign of g(x) is postulated to be the same for all values of x. 

iii. A and 6 are constants which measure the strength and range of scatterer. 
A is small in soft scatterers. 

In the case of a homogeneous spherical obstacle of radius 6 and propagation 
constant k h 




(r < 6), 

(r > b), 


The scattered wave function is given in (7) as a volume integral over the 
internal wave function. Hence if there is a scatterer whose g(x) approximates 
the required g(x) and whose internal wave function is known, then the internal 
wave function which corresponds to the approximate g{x) can be used in (7) 
to compute an approximate scattered wave function for the g(x) of interest. 

Let B-g a (r/a) be the function k 2 {r) - k\ which is associated with the approxi¬ 
mating scatterer. Then for a given scatterer of strength A and range 6 there 
is an approximating scatterer with fc 2 (r) - k 2 = B 2 g 0 (r/a) which is most 
closely related to it in the least square sense so that [14] 





'Poix)\ 2 dr = minimum. 


In the approximation method discussed in this and the next section we shall 
find scattered wave functions from (6) by using the internal wave function of an 
approximating scatterer whose parameters are chosen so that (21) is satisfied. 
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The simplest approximating function g 0 ( x ) that one might choose when 
\g{x)\ is nonincreasing is the step function: 


( 22 ) 


Then 


(23a) 

and 



9o(x) = 


A? 

0, 


AJ, 


L 

o, 


(if r < a), 
(if r > a). 

Or < 1), 

(x > 1), 


(236) 


£ 2 = 6] - ** 


> ^ if ^ ( b) > 0 ^ or x t 


^ 0 if g < 0 for all x. 

This approximation corresponds to computing the internal wave function by 

i eplacing the original scatterer by a homogeneous sphere. Since internal and 

scattered wave functions of a soft spherical scatterer are known [see equation 

( b)l one can substitute their analytical expression and (22) and (21) to 
obtain the values of k x and a which minimize (21). 

The wave function * 0 ( r ) which corresponds to B 2 g n (x/a ) [equation (22)] is 
[see equation (146)], to the first order in 6, — 6 0 , 

26 , 


^ exp (ik lT . so) 

= exp (ik 0 r • s,) + Ur), 


(H < a), 
(M > a), 

where Ur) is the scattered wave function. It is proportional to 6, - k 0 and 
for large r, varies as 1/r. Hence, for soft scattered we can neglect Ur) as 
compared with exp t6 0 (r • s„). Physically this corresponds to using the sphere 
internal wave function inside the approximating sphere and the Born approxi¬ 
mation for the “tail” of the scatterer, which lies outside the sphere. Hence 
we write 


2ft! 


(24) 


Mr) = kT+k 0 emU ''~ h, ex P (**> r -s.)> 


exp ik 0 (r • s 0 ), 


(l r l < a), 

(lr I > a). 


Equation (21), which defines a and k h becomes 


(25) / = 


46] 


(^1 + fto) 2 Jo 


(A! - A*) - A \f (L 


r 2 dr + A 4 



oo 


f 2 [^)r 2 clr 


I he equation 61/da - 0 and dl/dk u which determine the values of a and 6, 
which minimize (25), are 


(26a) 


AH 


= 26,(6] - 6«) 

Ski + A'o 
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and 


(266) 2fc?(fci + fro) 



<1 


0 


(fr! - frj) - A*g 



r 2 dr 


+ k 


' r. \ <* - 





r 2 dr = 0. 


In the limit as A*i — A 0 —► 0 and A 2 —> 0, the second integral in (266) becomes 
negligible, and (266) reduces to 


(27 a) 

where we define 
(276) 


-g(a)a :i = 



a 


x 2 g(x) dx, 


0 


a 

a ~ 6 


These equations relate the radius a of the approximating sphere to the range 
of the scatterer. In the limit, (26a) (which determines A'i) is 


(28) 


A'i - A'n — 2A' 2 f(a). 


6. Approximate Differential-scattering Cross Sections of Spherically Sym¬ 
metrical Scatterers. Once the parameters a and A*i of the approximating 
sphere have been obtained, the approximate scattered wave function for a 
spherically symmetrical scatterer (with the properties i, ii, and iii listed in the 
beginning of the last section) can be calculated by substituting \p(r) [equation 
(24)] into (6). 

The fundamental equation then becomes 


(29) 


UR) 


A 2 e lHk ° 

AtvR 



to(r)g ( g ) exp [ — ifc 0 (r •Si)] dr, 


where to(r) is defined by (24). After performing the angle integrations in the 
same manner as was done in the derivation of (166) and (176), we have 


(30a) + 9 (R) ^ A 2 R~'e iRk ° 


2A*i , 

- pi (k 1 

A i + A o / o 



oo 



r \ r* sin ra> . 
a i t 1 - dr 


+ 



r*g 


a 



r (x) 

sin ru 
ru 


— v 


sm ru 
ru 


(306) 

(30c) 

(30d) 


V = 


2 A* 2 


dr 


or 

u 


T , -r ex P - ko), 

M ~r /* o 

= A: 2 + A 2 — 2AiA 0 cos 0, 

= 21:o sin -£0. 


The approximate scattered wave function can be determined easily for a 
screened coulomb scatterer. In that case 


g(x) = x~ ] e~ x f 
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and after integration 
(31) \p 8 ~ A 2 e ikoR R~ l 


2 k x 


L*i + k 0 


0ta(ki—k o ) 


b 


b 2 e a/b sin u + bu cos ua 


1 + u 2 b 2 


The differential-scattering cross section is 


o) 1 -f- b 2 o) 2 

4 _ ^ p - a/i sin ua -f- bu cos ua 
u 1 + b 2 u 2 



(32) *.(B) = R*\+, | 2 . 

In the first approximation it can be shown [14] that 


a = | = 0.807, 

A 2 = ±0.379 (Ic 2 - A: 2 ). 

The corresponding results can be obtained with a Gaussian scatterer 


Then 

(33) * 
where 


$ 


(/Or) = exp ( -x 2 ). 

A 2 e iRk > j 2 ki bW ,, 7,3 

T I hTT, T+ 2 ‘-VW - >/(»)) 


j, x sin?/fl , ^ / 

^ ^ a - ^ 6+0 J 6 ^ C0S d x 

and the differential-scattering cross section is given by (32). Here in the first 
approximation [14] 

a = | = 1.235, 

±A 2 = 1.744(A 2 - k 2 ). 

The variables v, u, and « are defined by (30). As before, the differential- 
scattering cross section is given by (32). 

These expressions differ somewhat from those given by the authors in 
another paper [14]. In that paper it was assumed that 


MO = k~A~k l ex P M fc i - /co) + iki(r ■ so)], (for all r). 

The use of this approximation was based on the fact that f 0 (r) and <ty 0 (r )/dr are 
both continuous at r = a and /c 2 (r) — k 2 —> 0 rapidly as r — > oo. 

An error can be expected when the “tail” of the scatterer outside the 
approximating spheie extends out a large distance. One way of observing the 
existence of the error is to calculate the total-scattering cross section of the scat¬ 
ter as its range becomes large. In the case of the Gaussian scatterer the 
ratio of the total-scattering cross section to the “geometrical” cross section 



SCATTERING BY SOFT OBSTACLES 


115 


7r6 2 approaches zero as 6 —> » when one applies the formulas of the earlier 
paper [14]. This means that the tail is not given sufficient weight, for as long 
as a scatterer does not have a sharp boundary, this ratio should become infinite 
(but not as rapidly as it would in the Born approximation). We have given 
the tail more weight in the calculation made above by using the Born approxi¬ 
mation in the scattering by the tail. Although as the scatterer becomes large 
this overemphasizes the tail, formula (30a) should be valid over a considerably 
wider range than either the Born approximation or the simplified equivalent- 
sphere approximation used in the earlier paper [14]. 

Another approach to the problem of scattering by a soft obstacle has been 
made by Glauber [15]. He, independently of van de Hulst, has shown that 
equation (186) is valid for a wide variety of spherically symmetrical scatterers 
in addition to the homogeneous sphere. Glauber has derived (186) for the 
general case in several different ways. One is to replace the integral equation 
(6) by an approximate differential equation which is valid for very soft scat¬ 
terers. The internal wave function is found by solving the approximate 
differential equation. The scattered wave function is derived by substituting 
the approximate internal wave function into (6) and integrating. Also, by 
using the expansion of the scattered wave function in partial waves he improved 
(186) by replacing the term z = aOko by 2 = 2ak 0 sin 1-6. 

6. Total-scattering Cross Section for Very Soft Scatterers. The total¬ 
scattering cross section a can be obtained by integrating a s (6) over all angles, 

(34) a = 2 tt cr 8 (0) sin 6 dO. 

It is also related to the scattered wave function 

(34a) 

through the formula 

(35) 

This equation is equivalent to (34) only when (34a) is exact. In some approxi¬ 
mations (for example, in the Born approximation) the approximate f{6) is 
especially poor in the region of 6 = 0. Hence, even though it is much more 
convenient to use (35), the errors in the total-scattering cross section can better 
be identified with those of the/(0) as a function of 6 if (34) is used. Actually 
an alternative method of selecting parameters of approximating potentials is 
to make (34) and (35) and other similar formulas yield equivalent results. In 

the previous papers of this series total-scattering cross sections were derived 
from (34). 

Recently Glauber [15] developed a new formula for the total-scattering cross 
section of very soft scatterers which can be deduced without a detailed knowl¬ 
edge of the scattered wave function. It is a generalization of a result obtained 
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by van de Hulst [13] from (186) and (35). We shall give a derivation of his 
formula and discuss several special cases using it. 

Let us again apply the arguments of van de Hulst (see end of Sec. 3) to the 
scattering of a particle by a square-well (or step) potential. Our formula 
(146) for the internal wave function of a very soft sphere (as k—>0) shows 
that a particle which travels along a chord R is not deflected from its original 
path (as A’i » k 0 ) but that its wavelength is changed while it is in the range of 

the scatterer. The total change in phase 
by passage through the field depends on 
the distance u of closest approach of the 
particle to the scattering center. We let 
x(u) be this total change in phase and 
exp (izko) be the wave function of a free 
particle in the absence of the scatterer. 
Then at a large distance from the scatter¬ 
ing center, a particle which passes along a chord R has the wave function 

exp [izko + ixfa)]- 

The variation of the wave function due to scattering is 

gufro+«x(«) — f>izk 0 



so that the amplitude of the scattered wave is 



gizfro+i'xOO 


izk 


4 sin 2 £x(w). 


This result seems to be quite general. In the limit as the strength of a 
spherically symmetrical scatterer diminishes, there is no refraction but only a 
change in phase of the incident particle. Hence, in the limit the total-scat¬ 
tering cross section of a vanishingly weak spherically symmetrical scatterer 
becomes 



This expression, originally derived by Glauber, follows directly from the Glau- 
ber-van de Hulst formula by combining (35) with (186). 

The calculation of x(^) is easily performed when 



so that k(r) — k 0 
parameter is 


^(A 2 /k 0 )g(r/b) as A —» 0. In the limit as A —> 0, a useful 



(39) 
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We find x(bx) by noting that the change in phase per unit length at z along 
the chord /? is 


[A-(vV + z 2 ) - A-„] 


so that when (38) is satisfied 

x(u) = 



Hence c/irb 2 is given by (40) when 



(41) x(bx) 

and y is defined by (39). 

In some special cases 




/x- + v 2 ) dv, 


(a) g(x) = 



0 , 




Then 


x(bx) = 


2/(1 

0, 


(6) g{x) = e~ z ' 


~ x 2 )\ 

yields x(bx) = 


(c) g(x) = x 3 yields x(bx) = 


(d) g(x) = e 1 yields x(bx) = 

(e) g{x) = x~ l e~ z yields x(bx) = 


\y^e~ z \ 

V 
x 2 

yKi{x). 

2/A„(x). 


(x < 1), 
(x > 1), 


(if x < 1), 

(if x > 1). 


The function K,{x) is the jth Bessel function of the second kind of purely 
imaginary argument (see Watson (10, p. 78]). 

For a uniform spherical scatterer 



where H^y) is the Struve function [16] of order f defined by 


^ + 2 y 2 ) “ (sin y + y-' cos y) 

and y = 2ak 0 (m - 1) with m = ki/k 0 . This is exactly the result obtained by 

van de Hulst [13] and Hart and Montroll [12], As the range of thescatterer 
approaches co and hence as y —> oo j a/irb 2 —> 2. 
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In the case of the Gaussian scatterer 



where log 7 = 0.5772, and 

Ci x = log 72 ; — 

The function Ci x is tabulated in Janke-Emde. As x —» 00 , Ci x —> 0. Equa¬ 
tion (43) was first derived by Glauber [15]. 

ihe Born approximation of the total-scattering cross section is obtained 
from (40) by using the hypothesis that the phase change resulting from the 
passage of the incident wave through the scatterer is very small so that 
sin ix(bx) ~ ?x(bx). Then 





(1 - cos 0 n 







g(V(x 2 + v 2 ) dv 



Since y is proportional to b , the range of the scatterer, the Born approximation 
implies that as b —> <*, the total-scattering cross section increases with the 
fourth power of the range. This is incorrect for all scatterers considered. 
On the other hand, as b—> 0, the Born-approximation expression is in agree¬ 
ment with the results of more exact analysis. 

There is a case of a highly singular potential, g(x) = x~ z , for which the Born 
approximation (44) diverges and, even as b —■> 0, a does not approach 5 4 . 
In this case x(bx) = yx~ 2 so that a formal application of (40) yields 



There is some question as to whether equation (40) is valid for such a singular 
potential. Hence, equation (45) should be checked by a more rigorous analysis. 

The equation analogous to (40) for the total-scattering cross section per unit 
length of an infinite, soft, circularly symmetrical cylinder is 


(48) i - 2 sin- dz. 

As an example of an application of (46) let us consider the case of a uniform 
homogeneous circular cylinder of radius b. Then x(bx) is given by the case 
(a) under equation (41). Hence 

dx. 


2b 


= 4 



sin 


«a 

2 K 


x 2 )* 
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After letting x = sin 0, and integrating by parts, we have 

<r f x/2 

xr = 2y / sin 2 0 sin (y cos 6) dO = irHi(y), 

Jo 

where Hi(y) is the Struve function of order 1. This is a result previously 
obtained by van de Hulst [17] and Montroll and Hart [12]. 

7. Scattering by Systems Which Contain Several Scatterers. So far we 
have discussed only the scattering by isolated obstacles. In many physical 
applications of scattering techniques the influence of the scattering pattern 
of one scatterer on that of another becomes important. This influence can 
occur in several ways. The path length from source to observer may be dif¬ 
ferent in two of the scatterers (see Fig. 4a) so that a phase difference occurs 
between the two scattered waves and interference can take place. Also, 
multiple scattering can occur (see Fig. 46). The effect of multiple scattering 
increases with the concentration of scatterers. Here we shall be concerned 



only with the former effect. We shall assume that the concentration of scat¬ 
terers is so small that the multiple scattering is negligible. 

When the range of at least one of several different species in the same system 
is somewhat larger than that of the incident wave, the difference in the effec¬ 
tive path lengths of waves which pass through various species also contributes 
to the interference effect. This phenomenon is exaggerated as the difference 
in the average relative index of refraction of the various species increases. 
We shall now discuss the theory of this situation. The importance of the 
effect was first noted by Shomaker and Glauber [18] in their investigation 
of the theory of scattering of electrons by molecules containing heavy atoms. 

Consider a fixed distribution of atoms of a molecule in gas, liquid, or powder 
which are located at points {vj ] with respect to some preassigned origin. The 
scattered wave due to thejth particle is then, at a point r(|r| » |r ; |), given by 

(47) \pj = r~ l exp ( irk 0 )fj(6 ) exp [ifc 0 ( s 0 — Si) • r ; ]. 

Here s 0 and s> are, respectively, unit vectors in the direction of the incident 
and scattered radiation, cos 6 = s 0 *Si, and/>(0) is the function defined in (34a). 
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The sum of the individual scattered wave functions is the total scattered 
wave function 






exp [tfc n (s 0 - Si) • 


In gases, normal liquids, and powders the molecules are randomly oriented. 
Hence the total scattered intensity is obtained by averaging the intensity for 
a single molecule over all possible orientations. Then 


(49) 7(0) = r 2 (|^| 2 ) = V f* sm sr< \ 

jf ’’ No¬ 
where s = 2ko sin 4(9. 

r 'i = l r <- ~ Li is the distance between z'th and j'th atoms. 

When the scattering amplitudes are taken to be real (as they are in the Born 
approximation), the only interference effect is that due to the special distribu¬ 
tion of the atoms in a molecule. However, if the internal phase shifts differ 

foi different types of atoms, the interferences between waves scattered from 
these different atoms will be modified. 

the variation of the intensity of scattered radiation with angle in a single 
crystal is also obtained from (48). However, the simplification of random 
angular orientation which leads to (49) cannot be used. 

If we let 


(50a) 

then 

(505) 


fi = Pi exp i<p ]t 



1(B) = 



+ 



PjPk cos &jk 


j >k 


sin sr jk 
sr jk 


As a special application of this formula let us consider the scattering of 
electrons by a molecule of the form AB n , where the n B atoms are symmetri¬ 
cally spaced about the A atom. Let r n , r^, r " x , • • • represent the distances 
between the (B,B)-pairs and An, A' n , A'/n * • • the number of pairs at those 
distances. Also, let r 12 be the distance between an (A,B)-pair. Then 


(51) 1(6) - pi - 


np\ 


— 2pip 2 cos 8 12 


sin srj 2 
sr l2 


+ Pi l An + A' n + 


sr ii 


sr 


ii 


The combination cos 8 12 sin sr i2 might be written ^[sin (sri 2 — <$i 2 ) + sin (5r J2 
+ 5 12 )]. 

The procedure used to determine molecular structures by electron diffrac¬ 
tion is to guess a set of reasonable geometrical structures and sets of inter¬ 
atomic distances (the numbers r n , r' n , r[\, r n , • • • in the case at hand) and 
calculate a theoretical 1(6) for each of these models. That structure whose 
calculated 1(6) has maxima and minima at points closest to the observed ones 
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is considered to be the correct one. If the theoretical calculation is made 
using the Born approximation (as is usually done), the phase differences <5, ; 
are zero so that the locations of maxima and minima determined on that basis 
might not be the same as those given by (51). It can be shown that in molecules 
such as UF 6 which contain atoms of widely different atomic numbers, the 
phase shifts 8 {J may be sufficiently large even to change the sign of the first 
term on the right-hand side of equation (51). 

Let us choose the screened coulomb potential 


V(r) = Ze 2 r~ l exp (b = a 0 Z~^) 

as the scattering potential of an atom of atomic number Z. The constant a 0 
is the Bohr radius; a 0 = 0.528 X 10 -8 cm. Although the phase shift 6 i2 is 
angular-dependent, we can obtain its order of magnitude by considering only 
the first term inside the large parentheses of (31). In the notation of Sec. 4 
this potential is equivalent to 

yr 4 

g(x) = x~ 1 e z , A 2 = —2 me 2 — 

a o 

For an order of magnitude calculation of <5 i 2 we write 


612 ~ [a(k i - A*o )]a - [a(k x - /; 0 )] B , 

where the subscripts A and B refer to atomic species A and B. From the 
equations below (31) for the equivalent sphere radius a and index ki in terms 
of b and A 2 we can finally show (using the relativistic expression for k 0 ) that 




(X 0 /a 0 )(Z A — Z b) 
(E + £ 2 )i 


where X 0 /a 0 is the ratio of the Compton wavelength to the Bohr radius of the 
hydrogen atom. It is to be noted from (52) that as the energy of the incident 
electron increases, 5 i 2 —> 0, the Born-approximation value. As the difference 
in the atomic numbers increases, the phase <5 i 2 increases. Hence at very high 
energies and for atoms whose atomic numbers are close to each other, one is 
justified in using the Born approximation in electron-diffraction analysis. 

In Table I we have calculated approximate values of < 5 12 and cos 8 n for 
several molecules of the form AB n . It should be remembered that 5 i 2 is 
actually a function of angle, so that the values given are rough averages. For 


Table I 


Molecule 

6 i 2 

1 

cos 5 12 

CC1, 

0.4 

0.9 

SiF 4 

0.2 

1.0 

U F 6 

3.0 

-1.0 

WF 6 

2.5 

-0.8 

MoFc 

1 .3 

0.3 
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those values of cos 612 close to 1 the Born approximation is quite valid. On 
the other hand in the case of UF6, WF6, and M 0 F 6 , the peaks and minima are 
shifted considerably from the Born-approximation values, so much so that for 
many years the wrong structures were assigned to these molecules in order to 
make the theoretical and experimental scattering patterns fit. The impor¬ 
tance of the phase shift is the analysis of experimental electron-diffraction data, 

as was first pointed out by Shomaker and Glauber. The calculations in Table 

I were made for 40-kev electrons. 

0 

It can be shown that to a first approximation the effect of the internal phase 

shifts on crystals is to change the relative intensities but not the positions of 
the diffraction maxima. 

The authors are indebted to Dr. R. Glauber for several informative dis¬ 
cussions and for communicating some of his results to them before publication. 
Some of Glauber’s results have also been derived by Moliere [19]. 
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WAVE PROPAGATION IX HELICAL COMPRESSION SPRINGS 1 


BY 

E. H. LEE 

Abstract. In the past, spring surges have been analyzed on the assumption 
that adjacent coils do not collide during the motion. In many cases of com¬ 
pression springs subjected to impact this condition is violated. In this paper 
a theory of spring surges including coil closure is formulated, both for inelastic 
and elastic coil-impact conditions. Problems of the former type closely 
resemble certain problems in the propagation of plastic waves in compression, 
and the mathematical techniques developed for that work can be utilized. 
Perfectly elastic coil-on-coil impact demands a quite different type of mathe- 



Fig. l 

matical analysis associated with the propagation of force doublets along the 
elastic spring, resulting in intermittent pulses of compression. Some simple 
boundary-value problems are analyzed, with an indication of the generaliza¬ 
tion to more complex problems. The accurac} r to be expected from the type 
of analysis suggested is discussed. 

1. Introduction. We are concerned with the motion of a helical spring as 
depicted diagrammatically in Fig. 1. In practice such springs are commonly 
designed on the assumption of uniform force and strain distribution throughout 
the length of the spring, with a resulting linear relation between the transmitted 
force and the relative displacement of the ends. However, under conditions 
of dynamic loading, inertia forces of the spring itself modify the force trans¬ 
mitted, and wave-propagation effects occur. Such waves traveling through the 

1 The results presented in this paper were obtained in the course of research sponsored by 
the Ballistics Research Laboratory, Aberdeen Proving Ground, under Contract No. Da-19- 
020-ORD-798 with Brown University. 
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spring are known as spring surges, and they have been treated in the literature 
on the basis of the linear wave equation with constant velocity of propagation 
(see, for example, Dejuhasz [1]). The development of such a theory is briefly 
described below. In the case of helical springs loaded in compression, the 
applied force tends to close the coils together, and if adjacent coils come in 
contact, virtually no further compressive strain can take place at that section, 
khich closure will modify the motion based on the linear wave equation, and 
it commonly occurs in springs subjected to impact. The present paper 
describes a theory of spring surges taking into account coil closure. To our 
knowledge this is the first such treatment, previous work being limited to 
application of the linear wave equation to determine when coil closure first 
occurs [2], the subsequent motion not being analyzed. 

I he usual theory of spring surges is based on the consideration of motion 
of the spring wire parallel to the axis of the spring only, and the assumption 
that each element of the spring satisfies the force-longitudinal-strain relation 
of the whole spring. The complete analysis of the motion of a helix leads to 
the much more complex theory developed by Love [3], but experiment has 
shown satisfactory agreement with the simpler theory, as demonstrated, for 
example, by the work of Weibull [4]. We shall therefore base our development 
on the simple theory, but some consideration of the influences of the approxima¬ 
tions involved is given at the end of this section. 

In dev eloping the theory it is convenient to use as independent space variable 
the position x along the unstressed spring as depicted in Fig. 1. Such a 
Lagrange-type coordinate is used so that large displacements can be considered 
without complicating the wave equation which governs the motion. The 
displacement is u(x,t) 1 so that the position of an element of the deformed spring 
is given by x + u. Let / be the compressive force transmitted across the 
section x; then the equation of motion of an element dx is 

0 ) fz = —mu tl , 

where m is the mass of the spring per unit length and subscripts denote partial 
derivatives. Ihe nominal compressive strain e of an element of the spring 
is given by the reduction in length per unit initial length and is represented by 

e = -u x . Before closure occurs, there is a linear relation between the force 
and the strain, 

(2) / = E 6, 

where E is a constant depending on the dimensions and material of the spring. 
Substitution in (1) gives the linear wave equation for u, 

(3) u„ - ~ Uu = 0, 


with the constant velocity of wave propagation c 0 given by c 2 0 = E/m. Thus, 
as long as the linear relation (2) applies, spring surges travel with constant 
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(a) 


( 6 ) 


Fig. 2 


velocity when referred to the undisturbed spring coordinate x , even though 
the displacements may be large. The extensive theory of the linear wave 
equation can therefore be utilized. This type of analysis of spring surges is 
equivalent to the analysis of longitudinal waves in a linear elastic material. 
/ replaces the nominal stress, E the Young’s modulus, and m the density. 
Thus in applying this theory we in effect consider a spring as an equivalent 
elastic rod capable of withstanding large strains and having a low value of 
Young’s modulus and density, which depend both on the material of the spring 
wire and on the dimensions of the spring. 

In connection with the approximation mentioned above which is involved 
in this type of theory, that it is based on the longitudinal motion only and the 
static loading relation (2) for the whole 
spring, some details of the force transmission 
should be discussed. In loading a helical 
spring statically the line of force commonly 
acts along the spring axis so that a single p ** 
coil is loaded as shown in Fig. 2a. For Hat 
coiled springs, which include compression 
springs as discussed in this paper, it is suffi¬ 
cient to consider a single coil as a plane 
ring normal to the spring axis as depicted 
in Fig. 2, and the central loading of Fig. 

2 a then produces a constant torque PR in the spring wire and zero bending 
moment. The resulting constant twist produces a uniformly distributed dis¬ 
placement along the axis of the spring and thus a uniform strain. The inertia 
forces of the spring wire, however, act at the spring wire radius /?, so that their 
transmission through the spring must be similar to that of the force P shown in 
Fig. 2b. Such a force produces a varying torque and bending moment around 
each coil. If 6 measures the angular position of a section from the point of 
application of P, the torque isPP(l — cos d) and the bending moment PR sin d. 
These moment components will produce a different displacement pattern from 
that associated with axial loading. However, with inertia loading distributed 
around a coil, these moment variations will be averaged out. The torque 
PR(l — cos 0) has the average value PR and the bending moment PR sin 0 
the average value zero, so that the situation for axial loading is regained. 
Thus for a spring with many coils, if the average motion around a coil is of 
interest, and not the detail propagation of a stress wave within a single coil, 
the deformation associated with static axial loading can form a basis for the 
analysis. In many cases this averaging influence is enhanced by the circum¬ 
stance that the spring slides over a cylindrical former or inside a cylinder, 
which prevents bending of the spring axis and constrains the motion to be of 
the type associated with axial loading. Thus one would expect the simple 
theory to give satisfactory results, excluding the details of the deformation 
of an individual coil, and in fact it is found experimentally [4] to give good 



126 


E. H. LEE 


agreement even in the case of a wavefront of velocity discontinuity in which 
the averaging process is restricted by the rapid load gradients at the wavefront. 

The condition for coil closure is that the approach of sections of adjacent 
coils with the same angular position is equal to the initial separation of the 
coils. It is thus associated -with the difference in the displacements u at points 
separated in x by the pitch p of the spring coils and can be written 

( 4 ) u(x) — u{x + p) = p — d, 

wheie d is the axial dimension of the spring wire. Since our theory averages 
effects around a single coil, it is permissible to replace the finite difference by a 
derivative and so to consider coil closure to be associated with a fixed upper 
bound for the compressive strain e, determining the condition 

( 5 ) « = -«* < --- = (nun- 

V 

Such a maximum-strain condition forms the basis for the coil-closure analysis 
developed below. 

In general when closure occurs, there will be a relative velocity of impact 
between the adjacent coils, and the subsequent motion will depend on the 
rebound conditions at the point of contact. We shall adopt the usual theory 
of contact collision in which the relative velocity of approach and separation 
after bounce are governed by a coefficient of restitution. In this paper we 
shall deal with the two extreme cases: inelastic collision and perfect rebound. 
In the former case, which is dealt with in the next section, contact is maintained 
after collision until the surfaces are separated by the elasticity of the spring, 
there being no separating impulse generated by the collision. In the latter 
case of perfect rebound, the collision impulse which reduced the relative velocity 
of adjacent coils to zero is repeated to produce bounce with a relative velocity 

equal to that just prior to contact. As in the 
usual theory ( of collision, the relative approach 
during contact is neglected, for in general it will 
be small compared with the relative approach prior 
to closure. 

2. Closure with inelastic coil contact. As 

mentioned in the previous section, the force- 

strain relation for this case is that depicted in 

Fig. 3. The linear relation OA applies prior to 

closure, and during closure any increase of force 

T -> 0 max can be transmitted without change in strain. If 

1* ig. 3 & 

the force transmitted falls below the closure force, 
the linear relation OA in Fig. 3 again applies. This force-strain relation is a 
special case of a general functional relation 



( 6 ) 


/ = f(e) 
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which replaces the linear relation (2) of the discussion of the previous section. 
(3) then becomes 


(7) 


m 


Mix jf y>tt 0 , 


where/' is the derivative of /(e). Thus we have the wave equation in which the 
wave velocity c, given b}' c 2 = f'/m, is a function of the strain, i.e., of — u x . 
The treatment of boundary-value problems for this equation has been quite 
fully considered in the analysis of longitudinal plastic waves in rods, as devel¬ 
oped by Kdrmdn, Bohnenblust, and Hyers [5]. Thus, in effect, for the case 
of inelastic closure we can replace the spring by an equivalent nonlinear elastic 
rod of density m and with the stress-strain relation /(e) depicted in Fig. 3. 
The analysis developed in [5] is limited to cases in which the /(e) relation is 
concave downward, for which a wavefront of increasing stress spreads out as 



V 0 C D 

Fig. 4 


it propagates. The relation for a spring with closure depicted in Fig. 3 is 
concave upward, and the resulting modification of the theory, suggested by 
White and Griffis [0] and by Lee [7] must be adopted. In this case a wave of 
increasing stress may have a continually steepening wavefront, and a shock 
wave of force discontinuity may arise. 

Let us consider some simple boundary-value problems which provide 
examples of the types of behavior which this theory predicts. Figure 4 illus¬ 
trates the problem of a semi-infinite spring CD, 0 < x, initially at rest and 
unstressed, and subjected to a constant end velocity F 0 at zero time. If V 0 
is sufficiently small, as shown in Fig. 4, coil closure will not occur and the 
well-known linear elastic theory developed in the first section applies. A 
wave of constant material velocity V 0 and constant force and strain is propa¬ 
gated down the spring with velocity c 0 . The continuity equation determines 

(8) Vo = c 0 e 

and momentum considerations at the wavefront determine 

(9) / = mc 0 V o. 

This solution will apply if the strain given by (8) is less than the closure 
strain e mM . 
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For a larger value of prescribed end velocity, say V h (8) determines a strain 
greater than closure strain, and the solution shown in Fig. 5 applies. A wave- 
front of spring closure travels down the spring, producing conditions corre¬ 
sponding to the point B in the (/,e)-plane. Behind this front the coils are in 
contact and moving with the prescribed end velocity F,. The continuity 

equation and the known closure strain w determine the wavefront velocity 
c according to the relation 

( 10 ) Fi = <* m „, 

and momentum considerations at the wavefront determine the closure force by 

(11) / = mcV i. 

1 hus the velocity c of the closure wave increases linearly with Fi and the 
force linearly with (Fj) 2 . It is interesting to observe that the elastic part of 
the force-strain relation does not influence this solution as long as closure 




Vi C D 

Fig. 5 

occurs, but it does influence the magnitude of the impact at closure. The 
force required to overcome the elastic restraint in order to initiate closure 
transmits part of the impulse at the wavefront, and only the additional closure 
force is transmitted through the coil contacts. 

In the case of impact with velocity Vo, which does not cause closure, energy 
is conserved, and the work done at the point of impact is divided equally 
between kinetic energy and elastic-strain energy of the spring. In the case of 
impact causing coil closure depicted in Fig. 5, half the work done at the point 
of impact appears as kinetic energy; the rest represents work done on the spring 
and is equivalent, per unit length of spring, to the area OBG in Fig. 5. Only 
the portion OFG occurs as recoverable elastic energy, the part OBF being lost 
owing to the inelastic coil-on-coil collision. 

The present application of the theory of longitudinal waves in a rod with a 
concave upward stress-strain relation is more satisfactory physically than the 
previously available application of plastic waves in compression. The nominal- 
stress nominal-strain relation for a metal in compression is of the form shown 
in Fig. 6a. There is a vertical asymptote at the strain unity, since unit 
strain implies reduction of the length of the rod to zero, which would corre- 
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spond to squeezing the rod down to a flat sheet. The associated increase in 
area demands that the nominal stress increases without limit, so that a vertical 
asymptote occurs, and the stress-strain relation to he used in this type of 
theory has a concave-upward region. A longitudinal compressive wave in 
such a material, analyzed on the basis of longitudinal motion only, leads to a 
wave of stress and strain discontinuity, represented, for example, by the jump 
AB shown in Fig. 6a. Thus a wave of discontinuous change in section as 
shown in Fig. 66 is predicted. Such a wave would imply infinite lateral accel¬ 
eration due to the lateral expansion associated with the longitudinal compres- 





Fig. 6 




T~> 

IMG. / 

sion, and lateral-inertia effects would have an important influence in smoothing 
out such a wavefront, as indicated by the broken lines in Fig. 66, spreading it 
over several rod diameters. In the present application to spring-closure waves 
no such lateral-inertia effects arise, and this types of analysis can be expected 
to provide a closer approximation to the physical situation. 

The application to more complex boundary-value problems can be readily 
treated by considering the wavefront configurations in the (x,0-plane. Con¬ 
sider, for example, the problem shown in Fig. 7 of a constant end velocity V 0 
applied to a spring of length l, the other end of which is fixed. Suppose V 0 
does not cause closure because of the initial wave emanating from the impact 
end. Let A in Fig. 7 a represent the stress and strain values behind this 
initial wavefront OC in Fig. 76. Reflection at the fixed end at C based on a 
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linear elastic relation would produce a reflected wave of twice the stress ampli¬ 
tude of the initial wave, and the situation shown in Fig. 7 represents a case 
when such an elastic reflected wave would determine a strain greater than the 
closure strain. The reflected wave will then be a closure wave CD represented 

)y the j ump AB in Fi s* 7a - The conditions to be satisfied are that closure 
occurs behind the wavefront CD so that the particle velocity v there is zero to 

satisfy the fixed-end condition, and the continuity condition across CD becomes 

^0 = c (w - <u) = c 

Momentum considerations across the closure wave CD give 

fa — f a = mcV o = f B — mcoVo. 

(12) determines the wave velocity c, which gives the gradient of CD in Fig. 

76, and (13) then determines the closure force f B . When the closure wave 

reaches the impact end at D, the spring is closed completely and compatibility 

with the fixed-end condition demands that the velocity at the impact end 
should fall instantaneously to zero. 

More complex boundary-value problems can be handled most easily by 
using the characteristic relations in the (x, 0 -plane, 

( 14 ) f ± mcov = const, on x + c 0 t = const. 

when closure does not occur. Wavefronts of closure satisfy the continuity 
and momentum conditions, 

(15) Av = c At, 

Af = me Av, 

where A lepresents the discontinuity of the variable due to the passage of the 

wave and c is the wave velocity. Combination of these relations with the 

closure condition e = e max will determine the solution for prescribed boundary 
conditions. 

3. Closure with perfectly elastic coil on coil rebound. We shall now con¬ 
sider the case in which perfectly elastic rebound occurs if adjacent coils collide 
with a nonzero relative velocity of approach. The time of contact, which 
will be given by the Hertz theory of impact [ 8 ], will be extremely short com¬ 
pared with the time of wave propagation along the spring. Thus such impacts 
will occur as pairs of curves along one of which a moving-impact point force 
is acting and along the other an equal and opposite force. These curves will 
be separated by the pitch distance p in the ^-direction, since impact collision 
occurs in adjacent coils initially at a distance p apart. Since we are not con¬ 
sidering the detailed propagation around a single coil and we are replacing the 
actual coil-closure condition (4) based on the relative displacement of sections 
initially a pitch length apart by the maximum-strain condition ( 5 ), such pairs 
of forces can be replaced by a force doublet containing equal and opposite 
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forces moving down the spring in close proximity to each other. Apart from 
these force doublets, the spring motion will be governed by the linear wave 
equation as discussed in the first section, so that the analysis of a spring with 
bouncing coil closure is replaced by the analysis of a linear wave system with 
an appropriate distribution of external forces. In the (x,/)-plane, the solution 
will consist of the trajectories of force doublets, between which the linear wave 
equation will be satisfied. In order to visualize how such solutions are devel¬ 
oped, it is convenient to study an associated system with lumped masses. 

Let us consider a row of equal rigid balls of diameter d sliding freely in a 
straight groove and initially separated uniformly from each other by the dis¬ 
tance 5. Such a system can be considered as a simple analogue of a spring. 
The balls can approach each other through the distance <5, at which instant 


contact occurs and a force of any magnitude can be transmitted beween the 
balls. This model is equivalent to a spring with no elastic response before 
closure occurs. Let us consider a semi-infinite row of such balls and give the 
end one a constant velocity V leading to closure. For simplicity we first con¬ 
sider the case of inelastic impact discussed in the previous section. .After 
a time 6 V, collision between the first and second balls will occur, and there¬ 
after the second ball will continue to move with the velocity V, and the two 
balls will remain in contact. After a time n(d/V), where n is an integer, the 
(n + 1 )st ball will begin to move in contact with the first n balls. The wave- 
propagation velocity along the row of undisturbed balls is given by 



c — (d + 5) — 


since the wavefront progresses a distance d + 5 in each time interval 5/V. 
I his corresponds to the velocity based on the Lagrange coordinate x for the 
spring. In fact the closure strain is 



max 


6 

d + 5 


and (1(>) is a repetition of (10). In the present case of lumped masses, momen¬ 
tum considerations will determine a series of impulses to maintain the constant 
velocity of the first ball in place of the constant force for the case of contin¬ 
uously distributed mass. The motion of the system of balls is thus seen to 

be closely related to the analysis given in the previous section for inelastic 
coil collision. 

Let us now consider the case of the balls for perfectly elastic impact. If the 
first ball of the series is prescribed to move with the constant velocity V, after 
a time 8/V it strikes the second ball, which bounces forward with velocity 2V 
in order to maintain the same magnitude of the relative velocity on approach 
and separation. When this strikes the third ball after a further time interval 
5/2 F, the velocities of the two balls are exchanged, according to the behavior 
in collision of perfectly elastic balls of equal mass. Thus the second ball is 
brought to rest, and the third ball is propagated forward with velocity 2V. 
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This procedure is repeated, each ball being subjected to an impulse which 
projects it forward with velocity 2V and then a reverse impulse on striking 
the ball ahead which brings it to rest again. The first ball continues its for¬ 
ward motion with velocity V and again strikes the second ball, initiating a 
second double-impulse sequence which travels down the row of balls. The 
process is repeated as illustrated in Fig. 8, in which the abscissa is chosen in 



terms of the initial position .r and the displacement u to represent the balls as 
moving points. Since the velocity of the ball at the head of the wavefront is 
2V, the propagation velocity down the undisturbed row of balls is given by 


(18) c = 2(d + 8) 

O 

Behind this wavefront the balls are alternately at rest or in motion with the 
velocity 2F, apart from the first ball, which moves with the constant velocity 

V. Since the balls behind the wavefront are 
at rest and in motion during equal intervals 
of time, their average velocity is V and they 
move forward on the average in step with the 
first ball. Figure 8 shows clearly that the 
wave is propagated by a series of impulse 
doublets which move along the undisturbed 
line of balls with the velocity c of (18) and 
produce successively on each ball first a for¬ 
ward and then a backward impulse. This 
sequence provides a model for the considera¬ 
tion of the corresponding problem of spring surges. For the spring problem 
with a continuous distribution of mass, the moving impulse pair will be replaced 
by a moving force doublet, and we must investigate the response of the spring 
to such a force system. 

Figure 9 illustrates the stress or velocity distribution associated with a 
force doublet moving through an initially undisturbed spring. Since the 
doublet is produced by coil closure and rebound, it will move with a velocity 


max.) 
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c which is faster than the elastic wave velocity Co. Thus no disturbance can 
be propagated ahead of the forward force pulse EF in Fig. 9, and a wavefront 
of stress and velocity discontinuity arises. Since the forces are produced by 
closure impact, the strain in the region DE of coil contact must be the closure 
strain e max . Since we can replace the effect of coil closure by that of moving 
forces on a system which satisfies the linear wave equation, superposition can 
be used and the wavefront EF would simply be repeated with a negative sign 
to represent the motion due to the negative component of the force doublet 
at D . However, when we consider the associated force-strain relation for a 
spring with coil closure illustrated in Fig. 3, the concave-upward property 
suggests shock-wave conditions on loading and a spreading unloading wave 
having a tail propagated with the smaller velocity c n . The situation is there¬ 
fore examined in detail below, and it is shown that the region BC does in fact 
correspond to zero stress and velocity, so that nonzero stress and velocity occur 
only within the force doublet. Let f and v be the spring force and velocity 
along DE, and/i and v\ along BC. The momentum relation across the wave- 
front EF is similar to (15) but contains an additional term due to the external 
force T. It becomes 

(19) f+T = mcv. 

Continuity across EF is equivalent to the first of (15), 

(20) v = C€ max . 

But / is the elastic closure force mcle mhXf since additional contact forces are 
taken care of by the force doublet. Thus, for a prescribed value of v, T and c 
are determined. 

For the first part of the unloading wave DC, the momentum relation is 

(21) (/ - /j) + T = mc(o - vi), 

the continuity equation 

(22) (v — Vi) = c(c max - ei), 

and /1 and ei must be related by the force-strain law for the spring, 

(23) f ! = mc\t i. 

Subtracting ( 21 ) from (19) gives 

(24) f x = mcv i t 

and ( 22 ) from ( 20 ) 

(25) ^ = cc,. 

(24) and (25) combine to determine 

(26) /i = mc 2 d. 
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Thus, for c > c 0 , (23) and (26) are incompatible unless f x = e x = 0. Thus a 
force doublet produces a pulse of force, strain, and velocity in the spring which 
moves along with it, and B and C in Fig. 9 lie on the z-axis. Thus the rebound 
conditions on coil closure are such that the force and velocity conditions ahead 
of the closure are regained when it has passed. We note that as for the case 
of inelastic coil collision, the elastic force-strain relation does not influence the 
speed of propagation of the closure waves, but it does affect the magnitude of 
the impact force Z 7 . We shall now apply these conditions to simple boundary- 
value problems for springs with bouncing coils. 

We shall first consider a semi-infinite spring, x > 0, initially undisturbed, 
subjected to constant-velocity impact at the end x = 0. The impact velocity 



Fig. 10 Fig. 11 

V\ is sufficiently large to cause coil closure, so that substitution in (8) would 
indicate an elastic strain greater than e max . When the first coil of the spring, 
which is constrained to move with the impact velocity V x , strikes the next coil, 
the latter rebounds forward with velocity 2V\ because of the perfectly elastic 
assumption of coil-on-coil collision. This material velocity 2\\ determines 
the magnitude of the force doublet which represents the advancing closure 
front. The wavefront velocity c is given by (20) with v = 2V\ and the force- 
doublet magnitude by (19). This first doublet is followed by a sequence of 
identical force doublets generated by the constant-velocity boundary condition 
since the spring is left unstrained and at rest after each doublet has passed. 
The force and velocity distribution are therefore as depicted in Fig. 10, and the 
(a;,0-plane configuration is as shown in Fig. 11. Since the boundary velocity 
is V\ and the wavefront speed 2Fi/c max , the average strain behind the wave- 
front is e max /2. Thus the doublets must be spaced so that equal intervals of x 
are alternately strained to closure and unstrained. According to our theory in 
which the actual closure condition (4) is replaced by the maximum-strain 
condition (5), the strain pulses are infinitesimal in extent, depending on the 
Hertz time of impact. However, in the actual situation, the coil-collision 
forces are separated by a pitch length p in x , so that after an element of the 
spring has been set in motion by coil collision it will not be brought to rest 
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again until the wave has traveled a pitch length along the spring. This 
determines the scale of the velocity pulses illustrated in Fig. 10. 

Let us now consider a finite spring, 0 < x < l, subjected to constant-velocity 
impact at x = 0, and held fixed at x = l. Suppose'that the impact velocity 
V is such that closure first occurs on reflection at the fixed end. Figure 12 
shows the (x, 0-plane configuration. AO is the elastic wavefront, and the 
spring force and strain behind this front are represented by the point D in the 
force-strain diagram (Fig. 13). At A a reflected wave must be developed 
which will maintain the end of the spring at rest. The relative velocity of 
collision between the coil adjacent to A which is set in motion by the advancing 
elastic wave and the coil held fixed at A is V, so that elastic rebound determines 
a rebound velocity — V along the reflected wavefront AC. This magnitude 




will determine the force doublet due to closure and the speed of propagation 
of AC. With a relative velocity at the wavefront of 2F, continuity gives the 
wave speed c according to 

(27) 21' = c - 0 
Momentum considerations give the impact force according to 

(28) T +f E -f D = mc2V, 

where f E is the elastic closure force corresponding to E in Fig. 13 and f D the 
spring force in the region OAC of Fig. 12. The wavefront AC will be associ¬ 
ated with a doublet, behind which the conditions in OAC will be regained. A 
sequence of such doublets will be produced by reflection at the fixed end, and 
as in the previous case the average change in strain behind AC is half that at 
the doublet, so that regions of closure occupy half the spring length traversed, 
the intermediate intervals maintaining the conditions in OAC. Further reflec¬ 
tion will occur at the impact end, and the detail analysis of this will require 
a study of the interaction between the force doublets produced by reflection at 
each end, which is not treated in the present paper. More complex boundary- 
value problems can be treated by a development of this type of analysis. 
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4. Conclusions. The theory developed above forms a basis for the analysis 
of spring surges with coil closure for both inelastic and elastic collision. The 
two cases dealt with, inelastic and perfectly elastic collision, determine con- 
tlasting solutions, the rebound momentum due to coil bounce causes much 

higher wave-propagation velocities for the same impact velocity, the ratio 
being 2:1 for a simple closure wave in an initially undisturbed spring. While 
it is felt that the present type of analysis will be satisfactory, the assumptions 
involved in the consideration of simple longitudinal motion, which are dis¬ 
cussed in the introduction, cannot at present be assessed theoretically. Experi¬ 
mental check of the theory would therefore be very valuable. This is par¬ 
ticularly so in the case of elastic rebound, since the motion consists of pulses, 
each a coil pitch in wavelength, and it is on this scale that the detail of propa¬ 
gation round the helix may influence the motion. However, the fact that the 
simple elastic theory agrees with experiment [4] for a wavefront of velocity 
discontinuity, for which experimentally a wavefront exhibiting little dispersion 
was obtained, lends support to the suggestion that the type of theory developed 
in this paper will form a satisfactory method of analysis of spring surges with 
(‘oil closure. 
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ON THE WAVE EQUATION AND THE EQUATION OF 

EULER-POISSON 1 

BY 

ALEXANDER WEINSTEIN 

1. Introduction. A symposium on wave propagation is a suitable oppor¬ 
tunity to present some results on a remarkable class of hyperbolic partial dif¬ 
ferential equations which includes the classical wave equation as a special 
case. By considering this class of equations we shall obtain a new and per¬ 
haps the simplest solution of Cauchy’s problem for the wave equation. But 
it should be emphasized that this is not the sole purpose of presenting this 
paper. In fact, some of the equations considered here, such as the Tricomi 
equation, are known to have important applications in other fields such as 
transonic flow of compressible fluids. Moreover, our hyperbolic equations 
present a remarkable counterpart to a class of elliptic equations in generalized 
axially symmetric potential theory which also give rise to a number of applica¬ 
tions in hydrodynamics and the theory of elasticity. 

2. The Euler-Poisson-Darboux equation. Let x be an abbreviated notation 
for a system of m variables X\ } x 2 , * • • , x m , m > 1. In the following we 
write u(x,t) for u(x h x 2y • * * ,x m ,t). Let k be any real number, — oo < k < oo . 
We consider the hyperbolic differential equation 

( 1 ) Au k = u k tl + kt~ l u k , 

where the superscript k indicates the dependence of u on the parameter k . 
We write sometimes u (k) in place of u k . The symbol A is an abbreviation for 

m 

^ (d 2 /dx?). We intend to give here a solution of Cauchy’s problem for ( 1 ) 
1=1 

with the initial conditions 

(2) u k (x, 0 ) = f(x i,x 2 , • • • ,.r m ), wf(x, 0 ) = 0 , 

for a given sufficiently regular function /. The special case k = 0 yields a 
wave equation, A ?* 0 = Let us observe that, except for k = 0, equation ( 1 ) 
has variable coefficients. Moreover, we are considering here a singular initial- 
value problem since the coefficient kt~ ] is infinite on the hyperplane t = 0 . 
This implies that the usual existence and uniqueness theorems for Cauchy’s 
problems are not applicable to our case. 

For m = 1 , equation ( 1 ) appears already in Euler’s work. Actually, Euler 
considers an even more general equation containing, besides k, a second param¬ 
eter. Euler’s case was discussed in a famous paper by Poisson. A brilliant 
exposition of Poisson’s theory, always for m = 1 , was given by Darboux [ 1 ]. 
In our notation Darboux considered only the range 0 < k < 2, under the 

1 This paper was sponsored by the Office of Naval Research. 

137 



ALEXANDER WEINSTEIN 

mistaken impression that no essential new features would appear outside this 

range. Incidentally, the case k = £ corresponds to Tricomi’s equation. 

Besides the investigation of Euler’s equation Poisson has given in another 

memoir his famous solution of. Cauchy’s problem for the wave equation in 

ordinary three-dimensional space by using ( 1 ) for m = 3 and k = 2. It seems 

that the two cases m = 1 , 0 < k < 2 , and m = 3 , k = 2 , appeared unrelated 
in the classical work of Poisson. 

The cases m = 2,3, • • • , k = m - 1, and m = 1, ft = 1, 2, 3, • • • , were 
considered by Asgeirsson, who calls ( 1 ) the Darboux equation, a terminology 
which was letained in Courant-Hilbert [ 2 , p. 381], while Darboux, who in fact 
considers only the case m = 1 , calls ( 1 ) the Euler-Poisson equation. The cases 

m = k = _1 > “ 2 ; “3, • • • , were considered by M. H. Martin [3] and 
J. B. Diaz and M. H. Martin [4], in a remarkable extension of Riemann’s 
method. M. B. Kapilevic [5] has recently considered ( 1 ) for 0 < lc < 1 and 

m = 1 and m = 2 . Actually, the equation considered by Kapilevic, although 

seemingly more general, can be obtained from (1) by Hadamard’s method of 
descent [7]. The solution of Cauchy’s problem given by Kapilevic is very com¬ 
plicated and involves finite parts of integrals. In our exposition we shall 
consider Cauchy’s problem for all positive integers m and all real values of k 
and shall call ( 1 ) the equation of Euler-Poisson-Darboux (abbreviated EPD). 
The general method presented here is an improvement on that sketched in a 
preliminary note by the present author [ 6 ]. Further investigations were made 
by Diaz and Weinberger [9] and Blum. 

3. Fundamental recursion formulas. We have the following two funda¬ 
mental recursion formulas: 

(3) wffoO = tu k+ \x,t) } 

(4) u k (x,t) = t l ~ k u 2 ~ k (x,t). 

These recursion formulas allow us to obtain from a solution u k of equation ( 1 ) 
solutions of the same equation with the parameter k + 2 and 2 — k, respec¬ 
tively. Both formulas are mentioned in a somewhat different form in the 

book of Darboux, but only for the case m = 1 . As similar formulas occur in 

other parts of the theory of differential equations and as the proof is usually 
omitted, it is worth while to point out that such formulas are simple conse¬ 
quences of the following identities for functions v and w of a single variable t: 
Put w = t k ~ l v. Then 

(5) w tt + ~ ~ w t = t k ~ x {v lt + kt-h't). 

On the other hand, if we put tw = v t , we obtain 

, k + 2 

(6) w tt + —-— w t = 

In the following paragraphs we shall solve (1), (2) explicitly for k > m — 1. 
Then by a process which is essentially a repeated application of (3) and (4) 
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we shall obtain a formula for the solution of Cauchy's problem for the case 
k < m — 1 . 

4. Cauchy problem for the case k = m — 1 . It is well known [2] that the 

Cauchy problem for the EPD equation 

% 

(7) = i /"- 1 + m ~ 1 m 7 _1 

t 

with the initial conditions 



u m l (x,0) = f(x • • • ,x m ), = 0 


admits a unique solution given by a generalization of a celebrated formula of 
Poisson. 

Let 



= 2(x/w) m 

r(m/2) 


denote the surface of the unit sphere in the m-dimensional x-space. Then 
the solution of the Cauchy problem ( 1 ) and ( 2 ) is given by the formula 


(10) U rn l (xyt) - — j • • • J J(x 1 + OL\t,X 2 + Oi 2 t } • • • ,x m + OL m t) doj m 

= M(x,t;f), 

where M represents the mean value of / taken on a sphere of center x and of 
radius t. In (10) a l} a 2f * • • , a tn denote the components of a unit vector in 
x-space. The proof of our statement can be easily checked by direct computa¬ 
tions which are given explicitly in the book of Courant-Hilbert [ 2 , p. 411]. 
Sometimes we shall write in place of (10) in abbreviated notation 


u m ~' = M(x,t;f). 

6 . The method of descent and the cases k = m, m + 1 , • • • . Let us 

consider always for a fixed m the case in which the parameter k is one of the 

integers m, m + 1 , • • • . Let us imagine that for such a given k we have an 
equation of the type ( 1 ), viz., 



with the conditions 



( 12 ) u(x l} • • • ,x* + 1 | 0 ) = fi (xj, • • • ,x* +1 ), w # (x lf • • • ,x* + 1 , 0 ) = 0 , 

but in which we have, besides t, k + 1 independent variables xi, • • • , x k +i 

in place of Xi, • • • , x m . 1 his new Cauchyjjproblem is solved using ( 10 ) by 
the formula 


(13) u(xi, ■ ■ ■ ,x k+1 ,l) 


Wfc-J- 



j h(*i 


+ OC\t,X 2 + OC 2 t f • • • } Xk +1 + Otk+\t) d<jJk+l. 
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Let us put 


(14) 




j.t* m , 0. 


,0) = f(x 1 ,x 2 , 


Xm) i 


where / is the function which appears in the initial conditions (2). Let us 
moreover put 


(15) fi(x h x 2 , 


,x m , Xm+ij ' m * ,Xk+i) — fi(xi,X 2 , • * • , 3 :^, 0 , • • • , 0 ). 


In this way the u of (14) becomes a function only of x h x 2 , • • • , x m which 

satisfies equation (1) and the initial conditions (2) for the corresponding 
integral value of k > m. 

As the Cauchy problem (1) and (2) admits for such values of k only one 
solution (see Sec. 11), we can write in place of (13) 


(16) U = U k (x 1 ,3*2, * * • ,x m ,t) 

-Lf 

0>k +1 J 


J f(x 1 + ait, • • • ,X m + Omt) d(tik+l- 


The method described here for obtaining u k for k = m, m + 1 , • • • , is the 
celebrated method of descent introduced by Hadamard [7]. It is of interest 
for the following to rewrite the formula (16) by expressing the A>tuple surface 
integral as an w-tuple integral. For the convenience of the reader we shall 
sketch this reduction. As x and t appear in the fc-tuple integral of (16) only 
as parameters, it will suffice to consider the A*-tuple integral 



which obviously can be written by projecting on the equatorial hyperplane as 



f /(<*i, * • • ,a„,) da i • • • da m da m +i ‘ ‘ • da k 

J V(1 - ' - a£) - <4+1 

• • * — a Then (18) can be rewritten as 







»= i 


,a m ) da i • • • da m 



The second integral in (19) is of the same type as (18) for / = 1. Therefore, 
it is half of the area of a sphere of radius R in a space of k — m + 1 dimensions, 
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so that its value is ?R k m co*_ m +i. Putting this value in (19) we obtain 



COfc-m-f i 





With this result we obtain the following expression for (16): 


( 21 ) u k (x,t) = 


COA-f-l-m 



/(*r i + «i/, 


ori 2 + # - * +am 2 < l 


* £ m + Ot m t ) ( 1 CX ^ 



6 . The generalized method of descent. Formula (21) has been established 
for k = m, m + 1 , • • • . However, the integral on its right-hand side con¬ 
verges for all real values of k in the range k > m — 1 . We shall verify now 
that the function u k (x } t) defined by the right-hand side of ( 21 ) satisfies the 
differential equation ( 1 ) and the initial conditions ( 2 ) for all values of k which 
are greater than m — 1 . Of course, we have to substitute in (21) for the 
value 2(\/7r)yr(.$/2) as j n ( 9 ) even if s j s no t an integer. This extension of the 
solution to nonintegral values of k can be called the generalized method of 
descent. 

Let us put for a moment 


(22) U( x ,t) = u k (x,t). 

Then we have to consider in place of (21) the following formula written in a 
self-explanatory abbreviated form: 


(23) U(x,t) = f ■ ■ ■ I f(xi + <*,0(1 - p 2 ) (t - , -’" )/2 rfc* 1 • • ■ da m , 

(p 2 = a\ + • • • + a£). 

We have 



where /< denotes the derivative of / with respect to the argument x { + at. 
Integrating by parts, we obtain 




142 


ALEXANDER WEINSTEIN 


since each of the (m — l)-tuple integrals obtained in this process is obviously 
zero. Differentiating (25), we further obtain 


m 


(26) V 


7 “=1 


t 


m + 1 


i.j = 1 


Since 




+ Otit)<Xj{\ - p2) (h-m+l)/2 dai 


da m + 


Ut 

t 


^[«/0 - p 2 )(*-™+l)/2] = (1 - p2)«-m+l)/2 _ (k - m + 1) ctj (1 - p 2)(*—m—1)/2 

and since j = 1 , • • • , m, we obtain from (26) by an integration by parts 


m 


(27) + y r-i-rr 

t k — m 1 



i = i 



“I” 


[ —m(l - p 2 )(*-m+l )/2 + (k - m + l)a?(l - p 2 )(*-m-l)/ 2 ] ^ 


dotm> 


The m boundary integrals which should appear in (27) are easily checked to 
be zero. Using (25) and (27), we have 


m 


(28) U tt + 


'<'-11 

»=1 



filial + 


(! - p 2 )^-"* + 1 >/ 2 rf«i 


da m = AC/. 


Reverting back to as introduced by (21), we see from (28) and from the 
definition of U that u k indeed satisfies equations (1) and (2) for k > m — 1. 

Let us note that for (10) and (21) to be solutions it is sufficient to assume 
that / has continuous second derivatives. 

7. The case k < m — 1. In this section, it will be shown that, for k < 
in — 1 but k j* — 1, —3, —5, • • • , the Cauchy problem can be solved by 
means of the results of Sec. 6 and repeated application of the fundamental 
formulas (3) and (4). 

We begin by choosing a positive integer n such that k + 2n > m — 1. 
Using the methods of Sec. 6 , we solve the Cauchy problem for the function 
u k+ ‘ ln (x,t) satisfying the initial conditions 


(29) 


«* + 2 n (x\ 0 ) = 


/ 


j^m) 


(k + l)(k + 3) • • • {k + 2n - 1) 


> u^ 2n (x, 0 ) = 0 


By (4) we have 
(30) 

and, applying (3) n times, we obtain from (30) 


t k+2 


n— \yk+2n 


= U 


2—k—2n 


n 


(31) 


O_I- 

- 1/« A 


t dt 


( fk+2n-l u k+2n ) = u 
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Finally, applying (4) again, we obtain from (31) the fundamental formula for 
the solution of Cauchy’s problem, 





(J*+2"-l w *+2 ») 




It is clear from the above construction, which starts from u k+2n , that u k is a 
solution of the EPD equation. Further, by (32), 

(33) u k (x,t) = (k + 1 )(k + 3 )•••(* + 2 n - l)u k + 2 »(x,t) + Ctu k + 2n + 0(t 2 ), 

where C is a constant. Therefore, by (29), u k also satisfies the initial condi¬ 
tions (2). 

Let us recall that for u k+2n to be a solution it is sufficient that / have con¬ 
tinuous second derivatives. In order to be able to carry out the construction 
of u k , it is sufficient to require that / has n more derivatives. Thus, it is enough 
to assume that / has not less than \ (m — k + 3) continuous derivatives. 

If m — k is a positive odd integer, then n can be chosen so that 


k + 2n = tn — 1. 


Therefore, u k + 2n in (32) is given by a Poisson mean-value formula of the type 
(10). In this case, and only in this case, we have a Huygens' principle as already 
noted [9]. 

An important special case is the wave equation k = 0, for which (32) 
becomes a classical formula covering both even and odd ra’s. In both cases, 
for the wave equation it is sufficient to assume that / has not less than £(m + 3) 
derivatives. The assertion made [2, pp. 399/f.] that it is sufficient for / to 
have not less than £(m + 1) derivatives and its proof are incorrect. 

8. Representation of the solution by the Riemann-Liouville integral. 
According to a remark of Diaz and Weinberger, by introducing polar coordi¬ 
nates, solution (21), valid for k > m — 1, can be written in the form 


(34) u<H(x,t) 



- p 2 )(^ m— 1 )/2pm—1 dpj 


where M(x,t;f) is again a mean value as in (10). Incidentally, the analytic 
continuation of (34) as a function of k is the basis of the method used by Diaz 
and Weinberger to obtain a solution for k < m — 1 [9]. 

Let us now introduce the Riemann-Liouville integral 




(« - pY~ 1 F{p) dp, 


(35) 


H[F(s)} = 


o 


(0 > 0). 
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By making the substitution a- = pH 2 in (34) and then putting <* = 8 , we easily 
obtain the formula 



u {k) (x } \/~s) 





which is valid for k > m 1 . For k = m — 1 , (36) reduces to (10). 

On the other hand, as u k+u in (32) is given by a formula of the type (34), 
the same substitution allows us to write (32) in the form 


(37) u k (x,\/s) 



J (k+2n— m-f- 


l)/i [M (x,\/ s;/)s (m-2)/2 ], 


which is again valid for k > m — 1 — 2n, k ^ — 1 , —3, • • • 

9. Connections with the Asgeirsson mean-value theorem. Let us now 
assume that k = 0, 1 , 2, 3, • • • . We shall show first that the EPD equation 
can be considered as a Coulon, or ultrahyperbolic equation. 

Suppose for definiteness’s sake that k < m - 1. We consider a function 
u = n{x i, • • • ,x m ',y i, • • ■ ,y m ) which satisfies the Coulon equation 

(38) A z u = A„u. 

Then we have the following mean-value theorem [10]: 

(39) — / u(xi + ad,Q) du m = — f u(xi,ad) dco m . 

J com COm J com 

Putting u(x, 0) = f(x) and assuming that u depends only on x lt ■ • ■ , x m , y h 
■ • ■ , Vk+ i, we have by the method of descent (Sec. 5) 


(40) M(x,t;f) = 


./ m 



u(Xi,ait, 


• • 


.«*+i0 


X «* ^ 1 


» = i 


A+i 




(m-k- 3)/2 


a. 


da i 


da 


k+l 


%= 1 


If we further assume that u is radially symmetric in the space y h • • • , y k+h 

k+ 1 

and if we put p 2 = ^ a?, then (40) becomes 

i = 1 

(41) M{xi,t]f) = f u(xi, P t)(l - P *) d p. 

Um Jo 
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At the same time, equation (38) reduces to the EPD equation, 

A x u = u tt + ktrhi t , k = 0, 1, • • • , m — 2. 

Formula (41) has been considered for k = 0 by Asgeirsson. As in this special 
case, it is always an Abel integral equation the solution of which is given by 
formula (37). The difference between our consideration and Asgeirsson \s is 
that we assume that u in CouloiTs equation is radially symmetric in 
y u • • • , yi t+i, while Asgeirsson assumes that u depends only on y lt which he 
identifies with t. Needless to say, a similar procedure would lead to formula 
(30) for k = m — 1, m, m + 1, • • • . However, it should be realized that 
this procedure does not provide a proof for noninteger values of k and for 
negative values of/;. Even for k = 0, 1, 2, • • • , a verification of the formula 
would be required to make the application of Asgeirsson’s theorem independent 
of the methods developed in Sec. 7. Moreover, our method does not use Abel’s 
integral equations or analytic continuations. 

10. The exceptional cases of the odd negative integers. In the previous 
sections, the solution has been obtained for k 5^—1, —3, —5, • • • . The 
only assumption on / was that / should have continuous derivatives of order 
greater than or equal to ^(m + 3 — k). In that case, the solution also has 
continuous derivatives of the same order with respect to / even for t = 0. 
The situation is quite different in the exceptional cases A; = — 1, —3, —5, • • • . 
In these cases, a solution u (k) exists under the same differentiability conditions 
on /. However, the partial derivative d l ~ k u {k) / dt l ~ k of this solution becomes 
infinite like log / at / = 0 unless /(a*i, • • • ,.r,„) is polyharmonic of order 

id - k). 

In this generality, a solution was found by Diaz and Weinberger [9] and 
Blum in an unpublished paper, after the present author pointed out the excep¬ 
tional role played by polyharmonic initial values. Blum shows that a solution 
can be obtained directly from (32). B. Friedman was the first to give an 
example of a solution having nonpolyharmonic initial values, in a communica¬ 
tion to the present author. 

It is remarkable that for the exceptional values even a small deviation, no 
matter how regular, from polyharmonic initial values introduces logarithmic 
infinities in certain derivatives of the solution. In this sense, the solution 
does not depend continuously on the data. 

The exceptional role of polyharmonic initial values can be shown by the 
following elementary considerations: Let us first take k = —1 and assume 
that w^^O) exists. By putting / = 0 in 


= u\t 1) - 


u 


(-D 


t 


we find that Au (-1) (.r,0) = 0, which shows that f{x) must be harmonic 
For k = — 3, we have 


„(-3)( Xj 0) = lim ^ 

t-*o 


u ( ~ 3) 


t 
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From the EPD equation for k = -3, we see that 



Au (-3) (x,0) = lim 

*-♦0 



From (4) it follows that t l u\ 3) = ». If we require dV _3) /df 4 to exist 

for t = 0 and that the odd derivatives of u ( ~ 3) are zero, then d 2 u<~»/dt 2 also 
exists fort = 0. Therefore, Au<-»(x,0) = 0andby(42) wehaveAAw<- 3 >(x,0) = 0. 
This remark can be easily generalized to include all the exceptional values. 
In this case a solution of Cauchy’s problem is given by the formula 



u lk) (x,t ) = /(x) + 



/i=i 


__ A h f t 2h 

{k + l) ■ ■ ■ (lc + 2h - 1) 2 • 4 • • • 2h 


For arbitrary fix), let us give a brief sketch of Blum’s method. He observes 

that in virtue of (3) and (4) it is sufficient to consider the case k — 1. To 

show this, he first puts s = t 2 , as in Sec. 8, and gets in place of (3) and (4) the 
relations 


(44) u lk) (x,\/s) = « (fc+2) (x,\/i), 

(45) «<*>(*,\/s) = s"-»' 2 u 2 ~\x,y/s). 

Thus, if k = —(2r -f 1), repeated application of (44) and (45) yields 

(46) w-< 2r+1 > = As r+l u 2r+3 = As r+1 — 

ds r+1 

where A is a constant to be determined later. 

It is obvious that if the derivative of u a) in (46) is finite on s = 0, the result¬ 
ing solution u~ (2r+1) will be zero on s = 0. In fact, if (46) is to yield a solution 
w -( 2 f +i), w hich assumes nonzero initial values f(x h ■ ■ ■ ,x m ), a solution w (1 > 
should be constructed such that d r+l u^/ds r+1 has a singularity of the type 
s -(r+o q^is means that u (1) should have a singularity of the type log s. 
Such a solution is obtained by considering two appropriate solutions u w 
and u 2 ~ k for k = 1 + 2e and then passing to the limit as e approaches zero. 

11. Uniqueness. The uniqueness of the solution for k > 0 has been estab¬ 
lished by Asgeirsson by the method of Zaremba. However, for it < 0 there 
can be no uniqueness because any function of the type t l ~ k u 2 ~ k (x,t) which 
vanishes together with its derivatives with respect to t may be added to a 
solution u ik) (x,t) of the Cauchy problem. 

The situation has been recently clarified by E. K. Blum for the case m = 1 
in a paper presented to the American Mathematical Society. He showed 
that, for arbitrary k < 0, the difference of any two solutions is of the form 
t 1 - k u 2 - k (x,l), where u 2 ~ k {x,t) is a solution such that lim t~ k/2 u 2 ~ k {x,t) = 0 and 

lim t~ (k/2) ~ 1 u^~ k (x,t) = 0. 

(-*0 
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12. Conclusion. The key to the solution of all problems considered here 
was the identities (3) and (4) so long neglected. In conclusion, let us give 
another application of these identities. We define the Darboux operator, 
L a u = A u — v (t — ( a/t)u t . Then 

L a u & = uf = 

for all real a 1 s and 0’s. Therefore, for any values a iy • • • , a n we have 

Lf}+2nL an * * * L aj L ai u ( ^ 0. 

Taking /3 = 0 and «i = a 2 = • • * = a n = 2n, we obtain 

Lzi l u°(x i, • * * ,x m ,t) = 0. 

Putting m = 1, we obtain immediately a theorem of Friedrichs [2, pp. 416/f.]. 
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ON LlfiNARD’S DIFFERENTIAL EQUATION 


BY 


S. LEFSCHETZ 

The following differential equation 



x + 


dFjx) 

dx 


x + g(x) = 0 


has been investigated by Li£nard. He proved that under certain rather 
general conditions it possesses a unique periodic solution. A special case of (1) 
is the famous equation of van der Pol, which corresponds to F = ^(.t 3 /3 — x), 
g(x) = x. We propose to discuss an intermediary case, viz., the equation 



.. , dF(x) . , n 

x* + —-j— x* + x — 0. 
dx 


We shall assume that: (a) F{x) is an odd analytic function of x* for all real 
x’s and has only a finite number of zeros; (b) y = F(x ) behaves like a “ver¬ 
tical” parabola y = x n , n ^ 3; (c) from a certain finite positive value of x on, 
F(x) is monotone increasing with the obvious performance caused by (a) for 
x large negative. 

As usual, we introduce the associated pair of equations 

(3) x = y - Fix), y = -x , 

and conditions (a), (b), (c ) have for consequence (from Lienard’s treatment) 
that (3) has a finite number of limit cycles surrounding the origin. The 
outer limit cycle is stable outside. We have recently discussed the pattern 
of the solutions of van der Pol’s system outside the limit cycle (see a forth¬ 
coming paper in a monograph of the Annals of Mathematics Studies). We 
propose to extend the same results to the system (3). That is to say, we pro¬ 
pose to show that all the paths outside the outer limit cycle spiral clockwise 
around that limit cycle. 

Basically our study rested upon the behavior of the paths “at infinity.” 
Now the appropriate mechanism, which goes back in substance to Poincare, 
did consist in extending (3), or rather the equivalent equation 

(4) x dx + (y - F(x)) dy = 0, 

to the whole projective plane and examining the new critical points “at 
infinity.” That is to say, one replaces everywhere x by x/z, y by y/z , thus 
homogenizing the system. A point of the extended plane is now a triple 
(x,y,z) of numbers not all zero. Moreover ix,y,z) and (kx,ky,kz) represent 
the same point. The old points correspond to z ^ 0, and one may assume 
for them that z = 1; hence dz = 0. New points, however, are introduced, 
viz., on the line z = 0, and they correspond to the directions on the xy-plane. 
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Now for a system 

± = p ( x >y)> y = Q(x, y ), 

where P and Q are polynomials, I have shown that the extension is given by 



where n is the largest of the degrees of P and Q. The study of this relation 
along the line 2 = 0 proceeds practically as for the initial system. I have 
found that with a few minor modifications the same thing holds for the system 
(3). One verifies readily that the equation 






reduces to (4) upon making 2 = 1, dz = 0, and so it represents the appropriate 

extension of (4). The behavior of the paths of (3) at infinity is pretty largely 

governed by the behavior of the solutions of (5) about the critical points on 
2 = 0. 


Now it is found that, as in the case of van der Pol’s equation, the only 
critical points of (5) on 2 = 0 are point A at infinity on the .T-axis, or point 
(1,0,0), and point B at infinity on the y-axis, or point (0,1,0). 

Behavior at point A . We may take x = 1, dx = 0. The expanded equation 

becomes 



1 +y 





We are interested only in values of y,z near zero. This makes 



very large. Upon dividing by this expression, the terms which will have it in 
the denominator may be neglected, and what is left is, dropping the factor 2 , 
zdy = y dz. This integrates to y = Cz and shows that, as in the van der 
Pol case, A is the simplest type of node: the paths stream away from it in all 
directions. 

Behavior at point B. The full analysis of this point merely requires to 
prove the following property: there is just one path other than 2 = 0 tending 
to that point in the first quadrant. Now near B we may make y = 1, dy = 0. 
The equation thus becomes 



zx dx = 
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which ] 


M 



be written 



or again 



dx 

d log z 



F(x/z) 

x/z 



An examination of the solutions shows that they are parabolic, and in the 
xy-plane above the curve y = F(x). Let us suppose that we have two paths 
7 1; 72 tending to B and that 71 is above 72 . Then x/z—* + 00 along these 
paths, and F(x/z)/(x/z) is at least as large as (x/z) 2 . Let x x (z), x 2 (z) be the 
two solutions corresponding to 71 and 72 , and set x\ — x 2 = u > 0. Then 
from ( 6 ) 


du u f / n / m / x F(m) x 

(7) -T- } - = u - [<p(mi) - <p(m 2 ) ], <p(m) = ——, m = -• 

d log z X\X 2 , m z 

Now along 71 , y 2 almost 

s?(wi) = C ix pi , <p(m 2 ) = C 2 x pt , 

where either pi > p 2 or else pi = p 2 and C 1 > C 2 . Hence in (7) at the right 
the first term is small for 2 small, the second and third are negative, and the 
third is large. Hence for z small du/d log z < 0. Thus ?/ increases as z —> 0 
contrary to the assumption that 71 , 72 both —> B. Hence only one such path 
can exist. 

We have therefore proved our assertion that all paths outside the last limit 
cycle tend spirally to it. The signs of x,y show then that the spiraling is 
clockwise. 

We shall now examine another type of question closely related to the result 
just obtained. Consider the motion under a forced oscillation represented by 

/ox .. , dF(x) . 

(8) x H- x + x = e(t), 


where e is continuous and periodic with period r (hence bounded) and F(x) is 
as before. The question is to prove that there is a solution of ( 8 ) which is 
likewise periodic and of period r. The problem is evidently equivalent to 
finding a closed trajectory of 

(9) x = y - F(x ), y = -x + e(t) 

described in the time t. By a well-known procedure based on Brouwer’s 
fixed-point theorem, the problem reduces to finding a simple closed curve T 
at all points of which the vector IF, whose components are the right-hand sides 
of (9), points inward. Now it so happens that, in the literature, the generality 
of the system analogous to (9) causes V to be quite complicated. We propose 
to show that in the present instance T may be chosen exceptionally simple. 

Let generally |F|, • • • , denote the magnitude of the vector V, • • • . 
Let more particularly V be the vector whose components are the right-hand 
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sides of (3). 
trarily large. 


At points M(x,y) arbitrarily far from the 
For 


origin, |F| is arbi- 


l*T = (y - F{x)y + x 2 , 


and so either x is large, and \V\ likewise, or else x is bounded and y is large 
hence \V\ is again large. Since \ V — W\ = |e(*)| is bounded, to show that W 
points inward along T, it is sufficient to show that T is sufficiently far out and 
that V points inward (never tangentially) along r. 

To construct an appropriate r, we exploit an idea which the present author 

developed in his Princeton seminar of three years ago. It rests upon the 
consideration of the system 


x ~ (y F(x)) cos a + x sin a, 

y = (y ~ F(x)) sin a — X cos a, 

0 ^ a < tan a = m, 

whose vector V a is merely V rotated counterclockwise by the fixed angle a. 
The same method as before shows that (10) has perhaps a singular point at 
infinity at B (on the y-axis) and certainly at the point A a at infinity on the 
line y = mx. Since there are no asymptotes in directions in the second and 

fourth quadrants, all the solutions of (10) far out cross the x-axis. The curves 

(II) y = F(x) — mx , 

(III) y = F{x) + - 

m 

delimitate the regions x ^ 0, y ^ 0. Since far out they are in the first quad¬ 
rant, a solution of (10) starting at a point P on the x-axis far to the left behaves 
as in the figure and is arbitrarily far out, if P is remote enough. That is to 
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say, the arc PQ is arbitrarily far out if P is sufficiently remote. Now R is to 
the right of S where OS is the largest root of F(x) = rax, or the largest x of an 
intersection of y = F(x) with the line y = rax. Thus for ra sufficiently large, 
OS, hence OR, may be made arbitrarily large. Moreover increasing ra puts 
R even farther to the right. As a consequence, for ra and OP large enough 
the point R, and hence by symmetry the arc RT, will be arbitrarily far out. 
Hence the full solution of (10) from P to T will be arbitrarily far out. Further¬ 
more T is to the right of P, since otherwise the path PQRT extended below 
P would cross itself, which is ruled out. The arc PQRT of the solution of (10) 
plus the segment PT constitutes T. By construction along PQRT the vector 

V makes with the tangent to the arc directed forward a clockwise angle a 
and hence points inward throughout. Along the open segment PT the vector 

V points above the segment and hence again inward relative to F. Since T is 
arbitrarily far out, it satisfies all the required conditions, to establish that 
( 8 ) does have a periodic solution of period r. 

It may be observed that all the considerations of the present paper apply 
with only insignificant modifications if ( 1 ) is replaced, say by the more general 
Li 6 nard equation, 

.. dF(x) . , s n 

x H- :i — x + g(x) = 0, 


where g(x) behaves like F(x) and for x very large the order of g(x) is less than 
that of F(x). One will also replace ( 8 ) by 


x + 



x + g(x) = e(t), 


where e(t) is as before. 

Princeton University, 
Princeton, N.J. 




THE EFFECT OF SMALL CONSTRAINTS ON 

NATURAL VIBRATIONS 1 

BY 

R. J. DUFFIN AND A. SCHILD 

1. Introduction, lhe addition of a constraint to a vibrating linear mechan¬ 
ical system tends to increase the frequency of each normal mode. It is the aim 
of this paper to show how this qualitative statement can be expressed in a 
quantitative form. Consider constraints which can be varied continuously. 

We wish to find the small increase in frequency resulting from a small increase 
of constraint. 

lhe original system will be called the unconstrained system S; the system 
with the small additional constraints, the constrained system S'. Our problem 
is resolved by the following simple formula: 

(4) 6X = <57. 


Here X is the square of the angular frequency of a given normal mode of vibra¬ 
tion of S, and 6X = X' — X is the increase in X resulting from the increase in 
constraints. To maintain this increase in constraint, external forces of ampli¬ 
tude / must be applied to S ; the resulting vibration of S' is assumed to be 
suitably normalized. Now let the unconstrained system S be at rest. Let 
the forces / act to produce a static deflection. Then <57 is the potential energy 
of S in this state of static deflection. 

In many cases the normalized mode of vibration is known for the constrained 
system S'. The forces / can then be deduced. Hence <5X can be obtained 
without knowing any of the other constrained or unconstrained modes. 

lhe proof of formula (.4) can be carried out abstractly using the methods 
of Hilbert space. The essential hypothesis needed about S is that the kinetic- 
energy quadratic form 7 be completely continuous relative to the potential- 
energy quadratic form 1 . The complete proof will be given in another paper. 
Heie applications to plates, membranes, and other systems are sketched. 

Formula (.4) reduces the dynamic problem to a problem in statics. To 
facilitate the solution of this static problem, the following formula is employed: 


(B) 


8V = <57*. 


Let S* be a system which coincides with S in the neighborhood of the applied 
forces / of additional constraint. Outside of this neighborhood, S* is essen¬ 
tially arbitrary. Formula ( B) asserts that, under the same forces /, the poten¬ 
tial energies acquired by S and £* are equal. An equivalent way of stating 

1 Thls Paper is based on research conducted in part under a contract between the Flight 

° f thG U S ' Air F ° rCe and Carnegie In stitute of Technology, Contract 
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(B) is as follows: Close to a rigid constraint, the stiffness of an elastic body S 

is essentially independent of other constraints which are far away. Formula 

(B) is a general statement similar in nature to St. Venant's principle. In the 

applications S* is taken to be an infinite system because this simplifies the 
analysis. 

1 here are many applications of the basic problem considered. Various 
special cases have been treated in the literature by methods somewhat dif¬ 
ferent from those proposed here. It appears that these treatments lack unity 
and rigor. 

2. The clamped plate. The motivation for the present investigation was a 
nodal-line problem in the vibration of clamped plates. Consider a thin plate 
which is clamped on a closed curve B bounding a plane region R . If q is a 
deflection of the plate perpendicular to its surface, then in R the relation satis¬ 
fied b} r a normal mode of vibration is 



AAq = A q. 


Here A is the Laplacian operator, and the eigenvalue A is the square of the 
angular frequency of vibration. On the boundary B the clamping conditions 
are 



An important result in the theory of vibrations states that nodal lines cannot 
occur in the vibration of lowest frequency of a membrane. This was shown by 
Rayleigh. The question has been raised by Weinstein whether or not the 
same statement is true for the clamped plate. In other words, for the smallest 
value of A, is it possible that q will change sign in R? Weinstein's question is 
related to several other problems in vibrations. The one-dimensional analogue 
of the plate is the beam: it was shown by Rayleigh that the vibration of lowest 
frequency of a beam clamped at both ends cannot have a node. In a recent 
paper, Szego [3] has shown that, for clamped plates of a given area, the circular 
plate has the lowest frequency of vibration; he derives this result under the 
explicit assumption that nodal lines are not present. 

Weinstein's question has been investigated from three different lines of 
approach. In the first a semi-infinite strip is considered; its behavior suggests 
that there will be nodal lines in the lowest mode of vibration of a clamped plate 
in the shape of a long rectangle overlapping a circle of suitable radius [1]. In 
the second approach the doubly connected region between two concentric 
circles is considered; this eigenvalue problem can be solved explicitly and, for 
a sufficiently small inner circle, the lowest mode of vibration has a node along 
a diameter [2]. The third approach constitutes the main subject of this paper. 
This approach is more general and has applications to many problems of 
vibrations besides that of the plate. 
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Let R be a region, and let R be its mirror image in a tangent line l y which we 
choose as our x-axis (Fig. 1). Imagine a membrane or plate constrained on 
the boundary of R + R. Although there are actually two distinct regions, 
we can consider the eigenvalue problem for the whole region R + R simul¬ 
taneously. Corresponding to each eigenvalue X' there will be two eigenfunc¬ 
tions. One eigenfunction (<?'+) may be taken to be even with respect to the 
line /, the other (q'~) odd. Then 



q' + (z,y) = q' + (x,-y)> 

q'-(x,y) = -q'-(x,-y). 


Let us now open a small gap of length 2e between the two regions (Fig. 2). 
We shall call the eigenvalue problem with the gap open the unconstrained 




Fig. 1 


Fig. 2 


problem. Again the eigenfunctions will be even (</+) or odd ( q ~). But, in 
general, these will now have different eigenvalues. The constrained problem 
of Fig. 1 is degenerate, each eigenvalue X' corresponding to at least two eigen¬ 
functions. It is to be expected that opening the gap removes this degeneracy, 
so that each eigenvalue X' splits into two, X + corresponding to even states, 
X~ corresponding to odd states. This splitting of degenerate eigenvalues is, of 
course, a very common phenomenon. The Zeeman and Stark effects in quan¬ 
tum mechanics are well-known examples. 

The general formula ( A ) gives a first-order perturbation expression for the 
spectral fine structure. In the case of the clamped plate we shall see in Sec. 5 
that X~ is less than X + for small e. It follows that the mode of lowest frequency 
must be odd. The odd function has a nodal line at the gap. This is another 
counter-example which answers Weinstein’s question. 

3. General systems. We shall now give a short sketch of the method used 
to derive formula (A). The details of the proof will be published later. 

Consider the small vibrations of a general mechanical system S whose con¬ 
figuration is given by a vector q in a Hilbert space §. We can analyze such 
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a system in terms of a potential-energy quadratic form V(q,q) and a kinetic- 
energy quadratic form T(q,q). The amplitude q of a normal mode satisfies 

W V(q,x) - \T(q,x ) = 0 

for all vectors x in §. It is customary to normalize the eigenfunctions so that 

< 5 ) T(q,q) = 1 . 

Then it follows that 

(6) X = V(q,q). 

In fact, X may be defined as a stationary value of V{q,q) subject to (5). The 
stationary vector q satisfies the Eulerian equation ( 4 ). 

We now imagine that the system S is prevented from moving in a certain 
small region in the neighborhood of already existing constraints. This prob¬ 
lem with the small additional constraints, the constrained problem S', is 
characterized by equations similar to ( 4 ) and ( 5 ): 

(7) V(q',x') - \'T(q',x') = 0, 

(8) T(q',q') = 1. 

However, the vectors q', x' are restricted by the condition that they vanish in a 
subspace of using physical language, $ will be called the region of con¬ 
straint. More precisely, if is the orthogonal complement of SF in with 
the potential-energy form as metric, 

(9) F(S,£') = 0 , 

then we demand that q', x' be vectors in <£)'. In a certain sense, which we 
shall not specify here, the subspace must be small, corresponding to the 
physical fact that the additional constraints are small. 

For x in §, the expression V(q',x) — \'T(q',x) is a linear functional of x 
which vanishes in Therefore, there must exist a vector <f> such that 

(10) V(q',x) - \'T{q',x) = F( 0 ,x), 
and </> must be in tf, that is, 

(11) V(<t>,x') = 0. 

The configuration of the system described by </> is precisely that of the static 
deflection discussed in the Introduction. It is the static deflection of the 
unconstrained system S under forces / which are equal to the amplitudes of the 
additional forces of constraint acting in the vibration of S' when this vibration 
is normalized by ( 8 ). 

Subtracting (10) with x = q from (4) with x = q', we obtain 

(12) (X' - \)T(q,q') = -V(d>,q) 

= F(</>, 0 ) — V{(j),<t> + q). 
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For small constraints q and q f are close so that T(q,q') is essentially unity. 
Also, it can be shown that the term 7(0,0 + q) is of smaller order of magnitude 
than 7(0,0). Roughly speaking, this is due to the fact that in the region of 
additional constraints the kinetic energy is less important than the potential 
energy so that in the dynamic deflection q and the static deflection 0 are 
essentially equal in magnitude. Our sign conventions are such that 0 + q 
is then essentially zero in ^ and hence 7(0,0 + q) is small compared with 
7(0,0). Ignoring the small terms, (12) becomes 

(13) X' - X = 7(0,0), 

and this is, of course, our formula 

(A) 5X = 57. 

Formula ( A) holds in the simple form given here when the eigenvalue X' 
is nondegenerate. In the case of degenerate states the frequency changes are 
obtained by solving a secular equation. 

0 7T 0 7T~C 7T 

Fig. 3 Fig. 4 


Formula (A) is also valid with a dual interpretation of the potential energy 
57 in terms of the displacements of the original system in the region of the 
additional constraints. In this formulation, 5X can be obtained from the 
knowledge of a mode of vibration of the unconstrained s 3 r stem. 

4. The string. To illustrate our formula, we shall apply it to a simple 
problem where the solutions are known and where the approximate change in 
frequency given by (A) can be compared with the rigorous value. Consider 
a violin string stretched between points A at x = 0 and B at x = ir (Fig. 3). 
The violinist presses the string with his Anger at a distance e from one end 
(Fig. 4). This raises the tone. For a suitable choice of the physical units, 
the normalized nth mode of the constrained system AB' is given by 

( 14 ) g'n = 2 (tt - e)-i sill x, i f ({2 (lx = 1. 

7r — e Z Jo 

The force of constraint at B' is 



Under a force/„ at B' the static deflection 0 n of the unconstrained string AB 
is as shown in Fig. 5. The deflection at B' is 


( 10 ) 


< Pn(B f ) = €7T 1 (if ~ e)fn 
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The potential energy is therefore 

(17) SV n = MB')fn 

= 2671--*71 2 + 0(e 2 ). 

Formula (A) states that, correct to the lowest order in e, the change in the 
square of the angular frequency is 

(18) 5\ n = 2tw~ l n 2 . 

This result agrees with the known rigorous expression 

(19) 5X„ = n 2 ^1 - ^ - n 2 = 2tir-W + 0(e 2 ). 

To show how this calculation may be simplified by application of formula 
(B), consider the system S* in which the point A has been allowed to recede 


_ 

B' e B 

Fig. 6 

to x = — . Under the force f n at B' the static deflection <£* is as shown in 

Fig. 6, and 

(20) = ef n = 0„(B') + 0(e 2 ). 

Therefore 

(21) 8V* = bV + 0(e 2 ). 

The problem of the vibrating string is less trivial if the string has a variable 
distribution of mass along its length. The present methods still apply. 

5. The clamped plate (Continued). In the clamped-plate problem of Fig. 1 
we shall now make the simplifying assumption that the individual region R is 
sufficiently asymmetric so that its normal vibrations are nondegenerate. The 
problem of the combined region R + R is then twofold degenerate, but this 
degeneracy is due only to the symmetry of reflection in the line Z, and this 
symmetry is preserved when the gap 2e is opened. The even functions q + 
of the unconstrained problem tend to the even functions q' + of the constrained 
problem as e —> 0, and similarly for the odd functions. Therefore, the even 
and odd vibrations can be treated independently, and formula (.A) can be 
applied to each; this avoids a discussion of the secular equation mentioned 
near the end of Sec. 3. 

In the nth even vibration of the constrained plate the constraining force 
along the gap consists of a linear force density |dY„/dy 3 |. In the nth odd 
vibration the constraining force along the gap consists of a linear density 
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\d 2 (/ n /dy 2 \ of force moments. The functions <f n describe the deflection of the 
single region R when the gap is closed. Therefore, (f n and its derivatives are 
continuous along the gap and, since e is small, we can replace the force or 
moment densities of the constraining forces by their constant values at the 
center 0 of the gap: 



Let </>+ and <t>~ be the static deflections of the unconstrained plate under 
these respective force distributions along the gap; </>+ will be an even deflection, 
<t>- odd. Then the potential energy 5V n is given by 



Thus hV is determined by the behavior of the function <p in the immediate 
neighborhood of the gap. Formula ( B ) states that, to the lowest order in c, 
this potential energy will depend only on the part of the boundary at distances 
of order e from the gap and not on the shape of the boundary at finite distances 
from the gap. 

To make the argument more convincing, let us take the following example: 
Consider a large circular plate clamped around the rim. Here the Green’s 
function is known, and it can be verified that near the 
center of the plate the stiffness is low. This means that 
a small force acting there produces a large displacement; 
a small couple distributed along a unit segment near the 
center produces a large angular tilt. Let us now impose 
additional clamping along a diameter of the circle except 
for a finite gap at the center (Fig. 7). Now the plate is 
stiff in the neighborhood of the center: a small force act¬ 
ing there produces a small displacement, and a small 
couple distributed along a unit segment produces a small tilt. This shows that 
the stiffness near the gap is essentially determined by the clamping along the 
diameter and not by the boundary conditions on the distant rim. 

This physical discussion suggests that the principal part 8V* of the potential 
energy 8V can be determined from the static deflection of an infinite plate 
clamped along two half lines which are separated by a gap 2e as shown in 
Fig. 8. With this simplifying approximation the frequency increments can 
now be obtained explicitly, except for two constants which are independent 
of the geometry of our problem. Consider the standardized problem of an 
infinite plate clamped along two half lines which are separated by a gap of 
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length 2 as shown in Fig. 9. Let »+ be the potential energy under static 
deflection produced by a unit force density acting across the gap, and let tr 
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be the potential energy under a unit moment density across the gap. Then 
simple dimensional considerations show that 



IV? = et/>+, 

<5F? - = e 2 mlv~ 


Our formulas (-4) and (B) now state that, correct to the lowest order in 
the changes in the squares of the angular frequency are 



Subtracting these relations, we have 




Apart from the exceptional case when the second normal derivative d 2 (fjdy 2 
vanishes at the point of the boundary where the gap is located, the second 
term in (26) is negligible compared with the first for small values of e. Thus 
\ n is lower than X+. Physically our result shows that, for small gap lengths, 
the coupling between the regions R and R is due essentially to the rigidity of 
the plate, the gap acting as a fulcrum for a lever action. Pushing the region 
R down tends to cause the region R to rise. 

6. The membrane and acoustic vibrations. The membrane problem of 
Figs. 1 and 2 is similar to and simpler than the plate problem. For the mem¬ 
brane, = X' rigorously, and 



In this case it is easy to find an explicit expression for the constant v, since the 
static Green’s function for the region of Fig. 9 can be obtained by conformal 
mapping. 

Our formula (A) has an interesting acoustical application. Consider two 
rooms separated by a wall. If a small hole is cut in the wall, the frequency 
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changes can be computed. In this case one of the new frequencies is much 
lower than all of the original frequencies. This vibration can be regarded as a 
perturbation of a trivial zero-frequency mode of the two uncoupled rooms. 
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